This notebook contains all multivariate analyses of zoobenthic community structure using the new, nearly unheard-of modeling methods: packages mvabund, boral.
Again, to make it self-contained, there will be the same repetitive setup/data import/preparation part.
Setup!
Define the working subdirectories.
## print the working directory, just to be on the safe side
paste("You are here: ", getwd())
data.dir <- "data" ## input data files
functions.dir <- "R" ## functions & scripts
save.dir <- "output" ## clean data, output from models & more complex calculations
figures.dir <- "figs" ## plots & figures
Import libraries.
library(here) ## painless relative paths to subdurectories, etc.
library(tidyverse) ## data manipulation, cleaning, aggregation
library(viridis) ## smart & pretty colour schemes
library(mvabund) ## multivariate modeling analyses in ecology
library(boral) ## more multivariate modeling analyses in ecology
Organize some commonly-used ggplot2 modifications into a more convenient (and less repetitive) format. One day, I MUST figure out the proper way to set the theme..
## ggplot settings & things that I keep reusing
# ggplot_theme <- list(
# theme_bw(),
# theme(element_text(family = "Times"))
# )
## always use black-and-white theme
theme_set(theme_bw())
## helper to adjust ggplot text size & avoid repetitions
text_size <- function(text.x = NULL,
text.y = NULL,
title.x = NULL,
title.y = NULL,
legend.text = NULL,
legend.title = NULL,
strip.x = NULL,
strip.y = NULL) {
theme(axis.text.x = element_text(size = text.x),
axis.text.y = element_text(size = text.y),
axis.title.x = element_text(size = title.x),
axis.title.y = element_text(size = title.y),
legend.text = element_text(size = legend.text),
legend.title = element_text(size = legend.title),
strip.text.x = element_text(size = strip.x),
strip.text.y = element_text(size = strip.y)
)
}
## log y/min + 1 transform - useful for species counts/biomass data visualization
log_y_min <- function(y) {
log(y / min(y[y > 0]) + 1)
}
Import zoobenthic abundance data (cleaned and prepared).
zoo.abnd.sand <- read_csv(here(save.dir, "abnd_sand_orig_clean.csv"))
## convert station to factor (better safe than sorry later, when the stations are not plotted in the order I want them)
(zoo.abnd.sand <- zoo.abnd.sand %>%
mutate(station = factor(station, levels = c("Kraimorie", "Chukalya", "Akin", "Sozopol", "Agalina", "Paraskeva")))
)
Remove the all-0 species (= not present in the current dataset).
Maybe also remove the singletons (species appearing only once in the whole dataset and represented by a single individual = so rare that it’s unlikely they carry important information, but it would probably improve the run times).
(zoo.abnd.flt.sand <- zoo.abnd.sand %>%
select(-c(station:replicate)) %>%
select(which(colSums(.) > 0))
)
Perform a model-based unconstrained ordination by fiting a pure latent variable model (package boral - Hui et al., 2014). This will allow to visualize the multivariate stations x species data - similar to nMDS, can be interpreted in the same way.
I’m including a (fixed) row effect to account for differences in site total abundance - this way, the ordination is in terms of species composition.
NB this takes about a million years to run!
lvm.sand <- boral(y = zoo.abnd.flt.sand,
family = "negative.binomial",
## we want to control for site effects - there are 6 sites with 9 replicates each
row.eff = "fixed", row.ids = matrix(rep(1:6, each = 9), ncol = 1),
## 2 latent variables = 2 axes on which to represent the zoobenthic data
lv.control = list(num.lv = 2)
# ## example control structure, to check if function does what I want, because otherwise it takes an intolerably long time, and I'll shoot myself if I have to wait for it again
# mcmc.control = list(n.burnin = 10, n.iteration = 100,
# n.thin = 1)
#
)
Check the summary and diagnostic plots for the LVM.
summary(lvm.sand)
## model fit diagnostic plots
plot(lvm.sand)
The residuals plots look fine (no patterns in the residuals vs fitted, so variance is homogeneous, the quantile plot shows a normal distribution of the residuals) - the model fits the data pretty well.
Save the sand LVM.
write_rds(lvm.sand,
here(save.dir, "lvm_sand.RDS"))
Examine the biplot obtained by fitting the LVM, as well as the 20 most “important” species.
lvsplot(lvm.sand, jitter = T, biplot = TRUE, ind.spp = 20)
All in all, the final result resembles the nMDS ordination very much - same 4 clusters (Kraimorie + Chukalya, AKin, Agalina, Sozopol + Paraskeva). Kraimorie and Chukalya are better distinguished on the LVM plot than on the MDS, but still.
The run time is extremely, extremely long (~1h), but the data don’t need to be transformed, and the model fit can be examined and adjusted if necessary.
The species singled out as significant are probably somewhat different - have to check!
Redo the biplot, because this one is not very pretty. I’m not adding the species on top, first because I’m too lazy to figure out the procedure for ordering them, and second because the plot gets too busy.
## extract the LV coordinates of the stations from the model, so that the plot can be redone in ggplot
lvs.coord.sand <- as_tibble(lvm.sand$lv.median)
## add the stations from the original zoobenthic table (order was not modified)
(lvs.coord.sand <- lvs.coord.sand %>%
bind_cols(zoo.abnd.sand %>% select(station))
)
Make the plot and save it.
(plot.lvm.sand <- ggplot(lvs.coord.sand) +
geom_point(aes(x = lv1, y = lv2, colour = station)) +
scale_color_brewer(palette = "Set2", name = "station",
labels = paste0("S", as.numeric((unique(lvs.coord.sand %>% pull(station)))))) +
labs(x = "LV1", y = "LV2")
)
## save the LVM plot for the sand stations
ggsave(file = here(figures.dir, "lvm_sand.png"),
plot.lvm.sand,
width = 15, units = "cm", dpi = 300)
Let’s fit GLMs to the sites x species matrix to try and explain the observed differences in community structure by the variation of the environmental parameters.
These functions all come from package mvabund.
Import the environmental data - the one cleaned, prepared and saved in the previous notebook (classical multivariate methods). It contains long-term averages for the water column data (2009-2011 + 2013-2014) at each station, repeated for each replicate, and the sediment data (2013-2014), again repeated to the same number of replicates. Only the variables determined to be significant by PCA are kept.
Station is a factor, the rest of the variables are numeric.
Turn the zoobenthic data (minus the all-0 taxa) into a matrix - easier for the mvabund package and methods to deal with.
## there is already one subset of filtered count data (54 x 147) - use it
zoo.mvabnd.sand <- mvabund(zoo.abnd.flt.sand)
First, let’s see if the groups from the latent variable model (more or less equal to the clusters from the classical ordination) are valid, and which species exhibit a response.
## construct the vector of the clusters by hand, it's easier that way..
lvm.clusters.sand <- c(rep(1, times = 18), rep(2:4, each = 9), rep(3, times = 9))
## convert to factor
(lvm.clusters.sand <- factor(lvm.clusters.sand))
Check the model assumptions. 1. Mean-variance assumption => determines the choice of family parameter. Can be checked by plotting residuals vs fits: if little pattern - the chosen mean-variance assumption is plausible.
Another way: direct plotting (variance ~ mean), for each species within each factor level.
plot(manyglm(zoo.mvabnd.sand ~ lvm.clusters.sand, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.sand ~ lvm.clusters.sand, table = TRUE)
It’s not perfect, but it’s not too terrible either.
Everything looks more or less fine; fit the model.
glms.lvm.sand <- manyglm(zoo.mvabnd.sand ~ lvm.clusters.sand,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.sand)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.sand, which = 1:3)
### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
source(here(functions.dir, "default.plot.manyglm_grey.R"))
source(here(functions.dir, "plot.manyglm_grey.R"))
par(mfrow = c(2,2))
lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.sand, which = i, sub.caption = ""))
par(mfrow = c(1, 1))
I really don’t like the rainbow palette, but I would like to include these plots in my thesis results.. Will have to do something about it, just not right now.
Save the model!
write_rds(glms.lvm.sand,
here(save.dir, "glms_lvm_sand.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.sand.summary <- summary(glms.lvm.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
The factor (here - groups outlined by the LVM) is highly significant according to the models.
This also allows us to see which species exhibit a response to the chosen factor. The LR (likelihood ratio) statistic is used as a measure of the strength of individual taxon contributions to the observed patterns. I’ll save the summary for safekeeping, but I’ll also run an anova - to get an analysis of deviance table on the model fit (also better for extracting the species contributions, or at least I know how to do it).
write_rds(glms.lvm.sand.summary,
here(save.dir, "glms_lvm_sand_summary.RDS"))
Run the anova on the model.
(glms.lvm.sand.aov <- anova.manyglm(glms.lvm.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
I probably shouldn’t have printed all this out, but oh well who cares.
Save the ANOVA, too.
write_rds(glms.lvm.sand.aov,
here(save.dir, "glms_lvm_sand_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern (here - clusters in the LVM, which are really the different soft-bottom habitats).
top_n_sp_glm <- function(glms.aov, tot.dev.expl = 0.75) {
## helper retrieving the top n species with the highest contribution to the patterns tested by the GLMs, in decreasing order.
## Arguments: glms.aov - results from an ANOVA on the fitted GLMs
## dev.explained - proportion of explained deviance to use as cutoff
## get the change in deviance due to the tested pattern (= 2nd row from table of univariate test stats), and sort the species in order of decreasing contribution
uni.sorted <- sort(glms.aov$uni.test[2, ], decreasing = TRUE, index.return = FALSE)
## start at 10 species and check how much of the deviance is explained by their contributions. Repeat, increasing by increments of 10 until the desired explained deviance (set at function call) is reached.
top.n.sp <- 10
dev.expl <- sum(uni.sorted[1:top.n.sp])/sum(uni.sorted)
while(dev.expl < tot.dev.expl) {
top.n.sp <- top.n.sp + 10
dev.expl <- sum(uni.sorted[1:top.n.sp])/sum(uni.sorted)
}
## print the total deviance explained - just for information
print(paste("Total deviance explained:", round(dev.expl, 3)))
## return the final top species (and their univariate contributions, just in case)
top.sp <- uni.sorted[1:top.n.sp]
return(top.sp)
}
## get the top contributing species for the initial sand GLMs
(top.sp.glms.lvm.sand <- top_n_sp_glm(glms.lvm.sand.aov, tot.dev.expl = 0.75)
)
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.sand) <- names(top.sp.glms.lvm.sand) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.sand
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.sand <- zoo.abnd.sand %>%
select(station, names(top.sp.glms.lvm.sand)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.sand))))
)
plot_top_n <- function(top.n.sp.data, mapping, labs.legend, lab.y, palette) {
## helper for plotting top n species. Was hoping to avoid repeating it from way back when, but no dice.
## Arguments: top.n.sp.data - data frame (long) of top species' counts/biomasses at the different stations
## mapping - mappings of the aesthetics
## labs.legend - labels the use for the legend entries
## lab.y - custom label for y axis
## palette - custom colour palette (for consistency with other plots)
ggplot(top.n.sp.data, mapping) +
geom_point(alpha = 0.75) + # make points larger & partially transparent
scale_color_brewer(palette = palette, labels = labs.legend) +
ylab(lab.y) +
coord_flip()
}
(plot.top.sp.glms.lvm.sand <- plot_top_n(abnd.top.sp.glms.lvm.sand,
mapping = aes(x = species, y = log_y_min(count), colour = station),
labs.legend = paste0("S", as.numeric(unique(abnd.top.sp.glms.lvm.sand$station))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Well this is a nightmarish plot.. I’ll probably just put this awfulness in a table and call it a day, or play with lvsplot and the modeled ordination plot, if a plot is what’s needed.
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is hopelessly ugly and clumsy (and I’ll have to repeat it for the seagrass, too!), but I’m tired of being stuck on this. I still have many, MANY more things to do, and more time-consuming ones too..
top_sp_glms_table <- function(manyglms.obj.smry, group, p = 0.05) {
### extracts the top species in a group for which there is an observed effect in a manyglm test, at the specified probability level.
### Returns: tibble with the top species for the specified group/cluster, sorted (descending) by univariate LR value of the species, significant at the given p level.
## extract the univariate LR coefficients of the species and their p-values
sp_univar <- as_tibble(manyglms.obj.smry$uni.test, rownames = "species")
sp_p <- as_tibble(manyglms.obj.smry$uni.p, rownames = "species")
## combine in the same tibble
sp_all <- left_join(sp_univar, sp_p, by = "species")
## rename the columns
sp_all <- sp_all %>%
rename_at(vars(contains(".x")), list(~str_replace_all(., pattern = ".x", ".LR"))) %>%
rename_at(vars(contains(".y")), list(~str_replace_all(., pattern = ".y", ".p")))
## filter only the group/cluster we want, at the p-level we want
sp_all_flt <- sp_all %>%
select(species, contains(group)) %>%
filter_at(vars(contains(".p")), all_vars(. < p)) %>%
arrange_at(vars(contains(".LR")), list(~desc(.)))
}
top.sp.abnd.glms.lvm.sand <- lapply(names(glms.lvm.sand.summary$aliased), function(x) top_sp_glms_table(glms.lvm.sand.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.sand2" etc.
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% rename_at(vars(contains("lvm.clusters.sand")), list(~str_replace_all(., pattern = "lvm.clusters.sand", "group_"))))
top.sp.abnd.glms.lvm.sand <- lapply(top.sp.abnd.glms.lvm.sand, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.sand.long <- zoo.abnd.sand %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Kraimorie", "Chukalya") ~ 1,
station == "Akin" ~ 2,
station %in% c("Sozopol", "Paraskeva") ~ 3,
station == "Agalina" ~ 4))
## sum sp abundances by group; nest by group
zoo.abnd.sand.long.smry <- zoo.abnd.sand.long %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.sand <- map2(top.sp.abnd.glms.lvm.sand, zoo.abnd.sand.long.smry %>% pull(group), ~left_join(.x, zoo.abnd.sand.long.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.sand <- map2(top.sp.abnd.glms.lvm.sand, c(18, 9, 18, 9), function(x, y) x %>% mutate(mean_count = total_count/y))
)
To determine the relative taxon contribution to patterns: LR statistic - a measure of strength of individual taxon contributions. LR expresses how many times more likely the data are under one model than the other. This likelihood ratio, or equivalently its logarithm, can then be used to compute a p-value, or, compared to a critical value, to decide whether to reject the null model in favour of the alternative model.
In this case, the model shows which species exhibit a reaction based on the chosen groups - in other words, which species are more likely to be more/less abundant in each group.
For group 1 (= S1-S2), the species/taxa with significantly higher abundance are: Oligochaeta, H. filiformis, P. kefersteini, M. palmata, P. cirrifera, A. diadema (among others); and the ones with significantly lower abundance - even 0, in some cases - S. bidentata, B.lanceolatum, M. papillicornis, Melita palmata, P. jubatus, and so on.
For group 2 (= S3), the species with higher abundance are: B. lanceolatum, O. limacina, Oligochaeta (this is this strange artifact of 2013), P. kefersteini, L. flavocapitatus. The species with lower abundance are: H. filiformis, A. kagoshimensis, M. stammeri, Melinna palmata, etc. For group 3 (= S4-S6), the species with higher abundance are: C. gallina, L. mediterraneum - with very high dominance over practically all others; also Pseudocuma longicorne, Spio filicornis. The species with lower abundance are: H. filiformis, Oligochaetes (to a certain extent - they are still present, though), A. kagoshimensis, L. koreni, Harmothoe reticulata, Iphinoe tenella, Leiochone leiopygos.
For group 4 (= S5), the species with higher abundance are: Microdeutopus versiculatus, Eurydice dollfusi, Melita palmata, Polygordius neapolitanus, Polycirrus caliendrum, Polycirrus jubatus, Streptosyllis bidentata. The species with lower abundance are: A. kagoshimensis, Melinna palmata, P. cirrifera, P. ciliata, A. alba, I. tenella.
I love how the species with the highest variances (e.g. C. gallina, the most conspicuous example) are consistently pushed back - have lower LR scores. This is very good - C. gallina in particular is dominant in group 3, but is present also in all other groups - its substrate/depth preferences are very wide, so this is not uncommon. It’s not automatically pushed to the top of the list, but its reaction is detected by the manyGLM test. Neat! Contrast to the SIMPER results, where the species with the highest variance are consistently at the top - they contribute the most to the similarity, as per the test definition.
I’m going to save these as separate files (manually), then format them as tables - I know it’s a shame, but I’m too frustrated to figure out how to do it programmatically.
I’ll also put them in a word table in my final text, because I don’t want to deal with a million separate ones (embedded excel tables don’t split over multiple pages).
NB In my text, I’m switching the names/places of group 3 and 4, to be consistent with the SIMPER groups (I’m NOT going to repeat all this just to have the numbers match up). So the file names, table names, etc. remain as above. But in the text, I’ll have the following: group 1 = S1-S2, group 2 = S3, group 3 = S5, group 4 = S4-S6. REMEMBER THIS SO THERE IS NO CONFUSION!
Now, let’s try to see a different thing - which environmental parameters best describe the species response.
I’m going to use the PCA-filtered environmental data - it’s still going to be a slog, with 7 potential predictors..
First, construct the formula for the model - will do it separately in case I need to update it later, etc.
(formula.env.glms.sand <- formula(paste("zoo.mvabnd.sand ~",
paste(env.sand %>% select(-station) %>% names(), collapse = "+")))
)
zoo.mvabnd.sand ~ NH4 + NO3 + PO4 + seston + secchi + LUSI +
TOM + moisture_content + gravel + silt_clay
Fit the GLMs to the sand abundance data.
env.glms.sand <- manyglm(formula.env.glms.sand,
data = env.sand,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(env.glms.sand)
## all traditional (g)lm diagnostic plots
plot.manyglm(env.glms.sand, which = 1:3)
#
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.sand, which = i, sub.caption = ""))
# par(mfrow = c(1, 1))
Well, it’s good enough if you ask me (still the kinda strange “line” at lin.pred = -6; otherwise residuals are random enough).
Save the model!
write_rds(env.glms.sand,
here(save.dir, "glms_env_sand.RDS"))
Run the anova on the model - I want to see which predictors best explain the species abundance patterns I have. This is one function that would greatly benefit from being run in parallel..
(env.glms.sand.aov <- anova.manyglm(env.glms.sand,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 1.28 minutes...
Resampling run 200 finished. Time elapsed: 2.56 minutes...
Resampling run 300 finished. Time elapsed: 3.83 minutes...
Resampling run 400 finished. Time elapsed: 5.11 minutes...
Resampling run 500 finished. Time elapsed: 6.39 minutes...
Resampling run 600 finished. Time elapsed: 7.66 minutes...
Resampling run 700 finished. Time elapsed: 8.94 minutes...
Resampling run 800 finished. Time elapsed: 10.23 minutes...
Resampling run 900 finished. Time elapsed: 11.49 minutes...
Resampling begins for test 2.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 1.20 minutes...
Resampling run 200 finished. Time elapsed: 2.39 minutes...
Resampling run 300 finished. Time elapsed: 3.61 minutes...
Resampling run 400 finished. Time elapsed: 4.82 minutes...
Resampling run 500 finished. Time elapsed: 6.05 minutes...
Resampling run 600 finished. Time elapsed: 7.25 minutes...
Resampling run 700 finished. Time elapsed: 8.48 minutes...
Resampling run 800 finished. Time elapsed: 9.68 minutes...
Resampling run 900 finished. Time elapsed: 10.89 minutes...
Resampling begins for test 3.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 1.18 minutes...
Resampling run 200 finished. Time elapsed: 2.37 minutes...
Resampling run 300 finished. Time elapsed: 3.57 minutes...
Resampling run 400 finished. Time elapsed: 4.75 minutes...
Resampling run 500 finished. Time elapsed: 5.90 minutes...
Resampling run 600 finished. Time elapsed: 7.05 minutes...
Resampling run 700 finished. Time elapsed: 8.21 minutes...
Resampling run 800 finished. Time elapsed: 9.37 minutes...
Resampling run 900 finished. Time elapsed: 10.50 minutes...
Resampling begins for test 4.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 1.06 minutes...
Resampling run 200 finished. Time elapsed: 2.12 minutes...
Resampling run 300 finished. Time elapsed: 3.19 minutes...
Resampling run 400 finished. Time elapsed: 4.25 minutes...
Resampling run 500 finished. Time elapsed: 5.28 minutes...
Resampling run 600 finished. Time elapsed: 6.31 minutes...
Resampling run 700 finished. Time elapsed: 7.34 minutes...
Resampling run 800 finished. Time elapsed: 8.37 minutes...
Resampling run 900 finished. Time elapsed: 9.39 minutes...
Resampling begins for test 5.
Resampling run 0 finished. Time elapsed: 0.01 minutes...
Resampling run 100 finished. Time elapsed: 0.58 minutes...
Resampling run 200 finished. Time elapsed: 1.14 minutes...
Resampling run 300 finished. Time elapsed: 1.71 minutes...
Resampling run 400 finished. Time elapsed: 2.28 minutes...
Resampling run 500 finished. Time elapsed: 2.86 minutes...
Resampling run 600 finished. Time elapsed: 3.43 minutes...
Resampling run 700 finished. Time elapsed: 4.00 minutes...
Resampling run 800 finished. Time elapsed: 4.57 minutes...
Resampling run 900 finished. Time elapsed: 5.15 minutes...
Resampling begins for test 6.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.31 minutes...
Resampling run 200 finished. Time elapsed: 0.62 minutes...
Resampling run 300 finished. Time elapsed: 0.93 minutes...
Resampling run 400 finished. Time elapsed: 1.24 minutes...
Resampling run 500 finished. Time elapsed: 1.54 minutes...
Resampling run 600 finished. Time elapsed: 1.85 minutes...
Resampling run 700 finished. Time elapsed: 2.16 minutes...
Resampling run 800 finished. Time elapsed: 2.47 minutes...
Resampling run 900 finished. Time elapsed: 2.78 minutes...
Resampling begins for test 7.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.33 minutes...
Resampling run 200 finished. Time elapsed: 0.65 minutes...
Resampling run 300 finished. Time elapsed: 0.98 minutes...
Resampling run 400 finished. Time elapsed: 1.30 minutes...
Resampling run 500 finished. Time elapsed: 1.63 minutes...
Resampling run 600 finished. Time elapsed: 1.96 minutes...
Resampling run 700 finished. Time elapsed: 2.29 minutes...
Resampling run 800 finished. Time elapsed: 2.62 minutes...
Resampling run 900 finished. Time elapsed: 2.97 minutes...
Resampling begins for test 8.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.34 minutes...
Resampling run 200 finished. Time elapsed: 0.67 minutes...
Resampling run 300 finished. Time elapsed: 1.00 minutes...
Resampling run 400 finished. Time elapsed: 1.35 minutes...
Resampling run 500 finished. Time elapsed: 1.68 minutes...
Resampling run 600 finished. Time elapsed: 1.99 minutes...
Resampling run 700 finished. Time elapsed: 2.30 minutes...
Resampling run 800 finished. Time elapsed: 2.62 minutes...
Resampling run 900 finished. Time elapsed: 2.93 minutes...
Resampling begins for test 9.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.35 minutes...
Resampling run 200 finished. Time elapsed: 0.71 minutes...
Resampling run 300 finished. Time elapsed: 1.06 minutes...
Resampling run 400 finished. Time elapsed: 1.41 minutes...
Resampling run 500 finished. Time elapsed: 1.76 minutes...
Resampling run 600 finished. Time elapsed: 2.11 minutes...
Resampling run 700 finished. Time elapsed: 2.46 minutes...
Resampling run 800 finished. Time elapsed: 2.82 minutes...
Resampling run 900 finished. Time elapsed: 3.17 minutes...
Resampling begins for test 10.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.30 minutes...
Resampling run 200 finished. Time elapsed: 0.61 minutes...
Resampling run 300 finished. Time elapsed: 0.91 minutes...
Resampling run 400 finished. Time elapsed: 1.21 minutes...
Resampling run 500 finished. Time elapsed: 1.52 minutes...
Resampling run 600 finished. Time elapsed: 1.82 minutes...
Resampling run 700 finished. Time elapsed: 2.12 minutes...
Resampling run 800 finished. Time elapsed: 2.42 minutes...
Resampling run 900 finished. Time elapsed: 2.72 minutes...
Time elapsed: 1 hr 8 min 48 sec
Analysis of Deviance Table
Model: manyglm(formula = formula.env.glms.sand, family = "negative.binomial",
Model: data = env.sand)
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 53
NH4 52 1 1204.1 0.001 ***
NO3 51 1 481.4 0.001 ***
PO4 50 1 659.9 0.001 ***
seston 49 1 486.7 0.001 ***
secchi 48 1 913.4 0.001 ***
LUSI 47 1 -0.2 0.764
TOM 46 1 302.7 0.018 *
moisture_content 45 1 265.2 0.066 .
gravel 44 1 294.6 0.046 *
silt_clay 44 1 297.8 0.558
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Acanthocardia.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 2.178 0.999 1.738 1.000 3.445 0.987
Acari Actiniaria Alitta.succinea
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 7.262 0.449 1.738 1.000 1.214 1.000
Ampelisca.diadema Amphibalanus.improvisus Amphipoda
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 9.971 0.144 11.111 0.105 0.086
Amphitritides.gracilis Ampithoe.sp.
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.000 1.382 1.000 0.613 1.000
Anadara.kagoshimensis Aonides.paucibranchiata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 51.061 0.001 6.036 0.668
Apseudopsis.ostroumovi Aricidea.claudiae
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.313 1.000 7.654 0.375
Bathyporeia.guilliamsoniana Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.783 1.000 31.677 0.001
Bodotria.arenosa Brachynotus.sexdentatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.32 1.000 10.875 0.115
Brachystomia.scalaris Branchiostoma.lanceolatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.88 1.000 0.95 1.000
Caecum.armoricum Caecum.trachea Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 5.208 0.751 1.238 1.000 2.301
Capitella.minima Cardiidae Cerastoderma.edule
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.998 5.583 0.709 1.738 1.000 7.532
Cerastoderma.glaucum Chamelea.gallina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.399 0.714 1.000 40.253 0.001
Chondrochelia.savignyi Crassicorophium.crassicorne
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.656 1.000 0.046 1.000
Cumopsis.goodsir Cytharella.costulata Decapoda.larvae
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 3.445 0.987 16.488 0.010 4.145
Dinophilus.gyrociliatus Diogenes.pugilator
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.885 2.054 1.000 2.966 0.996
Donax.venustus Elaphognathia.bacescoi Eteone.flava
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 22.796 0.001 0.457 1.000 7
Eteone.sp. Eulalia.viridis Eunice.vittata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.526 2.055 1.000 0.058 1.000 2.984
Eurydice.dollfusi Eurydice.pontica
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.996 9.042 0.192 0.015 1.000
Exogone.naidina Gammaridae Gastrosaccus.sanctus
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.422 1.000 0.003 1.000 0.391
Genetyllis.tuberculata Glycera.sp.
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.000 4.617 0.827 1.83 1.000
Glycera.tridactyla Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 8.031 0.274 7.595 0.389
Harmothoe.reticulata Heteromastus.filiformis Hirudinea
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 12.862 0.038 60.52 0.001 4.652
Holothuroidea Hydrobia.acuta Hydrobia.sp.
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.825 3.138 0.995 3.582 0.971 1.932
Iphinoe.sp. Iphinoe.tenella
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.000 0.015 1.000 31.52 0.001
Kellia.suborbicularis Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 10.655 0.131 44.611 0.001 24.389
Lentidium.mediterraneum Leptosynapta.inhaerens
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.001 37.328 0.001 2.08
Lindrilus.flavocapitatus Loripes.orbiculatus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.999 0.996 1.000 10.096 0.142
Lucinella.divaricata Lysidice.ninetta Mactra.stultorum
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 27.256 0.001 0.015 1.000 4.196
Magelona.papillicornis Magelona.rosea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.878 32.322 0.001 6.909 0.570
Maldanidae Melinna.palmata Melita.palmata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 3.445 0.987 85.072 0.001 8.01 0.274
Microdeutopus.gryllotalpa Microdeutopus.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.009 1.000 3.445 0.987
Microdeutopus.versiculatus Micromaldane.ornithochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 6.59 0.595 16.373 0.010
Micronephthys.stammeri Microphthalmus.fragilis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.898 1.000 0.24 1.000
Microphthalmus.sp. Monocorophium.acherusicum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 6.038 0.668 8.355 0.240
Monocorophium.insidiosum Mysta.picta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 7.455 0.409 4.118 0.885
Mytilaster.lineatus Mytilus.galloprovincialis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 9.936 0.146 0.318 1.000
Neanthes.sp. Nemertea Nephtyidae
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.714 1.000 0.249 1.000 3.571 0.982
Nephtys.cirrosa Nereididae Nereis.pelagica
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.423 1.000 0.015 1.000 0.699 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 2.544 0.997 2.262 0.998
Nototropis.guttatus Odostomia.plicata Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.995 1.000 2.853 0.997 5.97
Ophelia.limacina Papillicardium.papillosum
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.673 0.002 1.000 3.477 0.987
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 15.816 0.012 1.351 1.000
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 21.34 0.003 5.282 0.739
Perioculodes.longimanus Pestarella.candida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 11.885 0.079 0.457 1.000
Pholoe.inornata Phoronida Phyllodoce.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 1.349 1.000 19.848 0.003 0.193 1.000
Pisces.larvae Pisione.remota Pitar.rudis
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 3.582 0.972 0.002 1.000 2.872 0.997
Platyhelminthes Platynereis.dumerilii
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 0.896 1.000 1.319 1.000
Polititapes.aureus Polychaeta.larvae
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 3.775 0.907 0.046 1.000
Polycirrus.caliendrum Polycirrus.jubatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 4.988 0.783 7.767 0.341
Polydora.ciliata Polygordius.neapolitanus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 42.608 0.001 2.208 0.999
Prionospio.cirrifera Protodorvillea.kefersteini
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
NH4 47.902 0.001 0.422 1.000
Protodrilus.sp. Pseudocuma.longicorne Rapana.venosa
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
NH4 0.059 1.000 22.859 0.001 16.348
Retusa.truncatula Retusa.variabilis
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
NH4 0.011 1.738 1.000 6.909 0.570
Rhithropanopeus.harrisii Rissoa.membranacea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
NH4 0.457 1.000 5.749 0.700
Sabellaria.taurica Salvatoria.clavata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
NH4 16.348 0.011 2.667 0.997
Schistomeringos.rudolphi Sphaerosyllis.hystrix
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
NH4 5.67 0.700 1.732 1.000
Spio.filicornis Spionidae Spisula.subtruncata
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
NH4 13.886 0.026 0.322 1.000 38.478 0.001
Stenothoe.monoculoides Sternaspis.scutata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
NH4 0.457 1.000 0.059 1.000
Steromphala.divaricata Streptosyllis.bidentata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA>
NH4 1.758 1.000 7.782 0.339
Syllis.gracilis Syllis.hyalina Tellina.fabula
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA> <NA>
NH4 0.088 1.000 5.847 0.688 7.655 0.375
Tellina.tenuis Tritia.neritea Tritia.reticulata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept) <NA> <NA> <NA> <NA> <NA>
NH4 4.644 0.827 21.326 0.003 4.162
Turbellaria Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept) <NA> <NA> <NA> <NA> <NA>
NH4 0.885 0.639 1.000 0.457 1.000
[ reached getOption("max.print") -- omitted 9 rows ]
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
The results suggest that the long-term average water column parameters have a major influence on the observed community structure; also the sediment TOM, and (marginally) the sediment composition (gravel content).
Save the ANOVA - I really, really don’t want to have to repeat it.
write_rds(env.glms.sand.aov,
here(save.dir, "glms_env_sand_anova.RDS"))
Get the taxa with the highest contributions to the tested pattern (here - species most affected by changes in water/environmental quality parameters).
(top.sp.glms.env.sand <- top_n_sp_glm(env.glms.sand.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.783"
Melinna.palmata Heteromastus.filiformis Anadara.kagoshimensis
85.071520 60.520034 51.060684
Prionospio.cirrifera Lagis.koreni Polydora.ciliata
47.902260 44.610886 42.608041
Chamelea.gallina Spisula.subtruncata Lentidium.mediterraneum
40.252681 38.478437 37.328312
Magelona.papillicornis Bittium.reticulatum Iphinoe.tenella
32.322228 31.676928 31.519717
Lucinella.divaricata Leiochone.leiopygos Pseudocuma.longicorne
27.255927 24.389432 22.859494
Donax.venustus Parvicardium.exiguum Tritia.neritea
22.796122 21.339562 21.326460
Phoronida Cytharella.costulata Micromaldane.ornithochaeta
19.847505 16.488442 16.372677
Rapana.venosa Sabellaria.taurica Paradoneis.harpagonea
16.348207 16.348207 15.816268
Spio.filicornis Harmothoe.reticulata Perioculodes.longimanus
13.885810 12.861565 11.884939
Amphibalanus.improvisus Brachynotus.sexdentatus Kellia.suborbicularis
11.110812 10.875153 10.654626
Loripes.orbiculatus Ampelisca.diadema Mytilaster.lineatus
10.096263 9.971170 9.935577
Eurydice.dollfusi Monocorophium.acherusicum Glycera.tridactyla
9.042071 8.354629 8.031053
Melita.palmata Streptosyllis.bidentata Polycirrus.jubatus
8.010391 7.782071 7.767249
Tellina.fabula
7.655439
I’m going to plot these top contributing species, but I’m not using the plot. At least this time it’s more manageable, but still not presentable enough..
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.env.sand <- zoo.abnd.sand %>%
select(station, names(top.sp.glms.env.sand)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.env.sand)))) %>%
## add clusters from LVM as a column
mutate(group = case_when(station %in% c("Kraimorie", "Chukalya") ~ 1,
station == "Akin" ~ 2,
station %in% c("Sozopol", "Paraskeva") ~ 3,
station == "Agalina" ~ 4))
)
(plot.top.sp.glms.env.sand <- plot_top_n(abnd.top.sp.glms.env.sand,
mapping = aes(x = species, y = log_y_min(count), colour = factor(group)),
group = abnd.top.sp.glms.env.sand %>% pull(group),
labs.legend = unique(abnd.top.sp.glms.env.sand$group),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Maybe a bit later I’ll try to get this nightmare above as a table…
Final analysis to try: which species respond differently to different environmental parameters? (= traits analysis - fit single predictive model for all species at all sites, but w/o attempting to explain the different responses using traits - the species ID is used in place of a traits matrix).
NB only use the top species that exhibited a reaction in the environmental model fit (= the ones accounting for ~75% of the total variability), and only the significant predictors - to improve run times.
When using LASSO (method = “glm1path”), the algorithm fails to converge - I’m not sure how to interpret it.. Maybe because the function tests each individual species:env.parameter interaction (does it really??), and none of them by themselves are sufficient to explain a species’ response. Not to mention the fact that the samples are not really independent (they are replicates at 6 sites, repeated 3 times).
When using method = “manyglm”, the result is the one shown above. It’s still a bitch to interpret - for example, what is the interpretation of an increase in abundance with both high PO4 and Secchi? Or with high NH4, but low NO3? Where are these conditions ever met?
In fact, everything points towards the conclusion that a species response is determined by a combination of eutrophication parameters in its environment (water column characteristics), and the composition of the sediments (organic matter and granulometry).
This is actually exactly the same thing that the PERMANOVA gives, in this particular case. However, in the future, I’m leaning more towards the modeling approach - it allows you to check the model fit to one’s real data; also, there are no data reductions due to calculation of distance matrices.
Import zoobenthic abundance data (cleaned and prepared).
zoo.abnd.zostera <- read_csv(here(save.dir, "abnd_zostera_orig_clean.csv"))
Parsed with column specification:
cols(
.default = col_double(),
station = [31mcol_character()[39m,
habitat = [31mcol_character()[39m,
replicate = [31mcol_character()[39m
)
See spec(...) for full column specifications.
## convert station to factor (better safe than sorry later, when the stations are not plotted in the order I want them)
(zoo.abnd.zostera <- zoo.abnd.zostera %>%
mutate(station = factor(station, levels = c("Poda", "Otmanli", "Vromos", "Gradina", "Ropotamo")))
)
Remove the all-0 species (= not present in the current dataset).
Maybe also remove the singletons (species appearing only once in the whole dataset and represented by a single individual = so rare that it’s unlikely they carry important information, but it would probably improve the run times).
(zoo.abnd.flt.zostera <- zoo.abnd.zostera %>%
select(-c(station:replicate)) %>%
select(which(colSums(.) > 0))
)
Perform a model-based unconstrained ordination by fiting a pure latent variable model (package boral - Hui et al., 2014). This will allow to visualize the multivariate stations x species data - similar to nMDS, can be interpreted in the same way.
I’m including a (fixed) row effect to account for differences in site total abundance - this way, the ordination is in terms of species composition.
NB this takes about a million years to run!
lvm.zostera <- boral(y = zoo.abnd.flt.zostera,
family = "negative.binomial",
## we want to control for site effects - there are 6 sites with 9 replicates each
row.eff = "fixed", row.ids = matrix(rep(1:5, times = c(8, 8, 4, 8, 4)), ncol = 1),
## 2 latent variables = 2 axes on which to represent the zoobenthic data
lv.control = list(num.lv = 2)
# ## example control structure, to check if function does what I want, because otherwise it takes an intolerably long time, and I'll shoot myself if I have to wait for it again
# mcmc.control = list(n.burnin = 10, n.iteration = 100,
# n.thin = 1)
#
)
module glm loaded
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 3008
Unobserved stochastic nodes: 3452
Total graph size: 21760
Initializing model
|
| | 0%
|
|++++ | 8%
|
|++++++++ | 16%
|
|++++++++++++ | 24%
|
|++++++++++++++++ | 32%
|
|++++++++++++++++++++ | 40%
|
|++++++++++++++++++++++++ | 48%
|
|++++++++++++++++++++++++++++ | 56%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|* | 3%
|
|*** | 5%
|
|**** | 8%
|
|***** | 11%
|
|******* | 13%
|
|******** | 16%
|
|********* | 19%
|
|*********** | 21%
|
|************ | 24%
|
|************* | 27%
|
|*************** | 29%
|
|**************** | 32%
|
|***************** | 35%
|
|******************* | 37%
|
|******************** | 40%
|
|********************* | 43%
|
|*********************** | 45%
|
|************************ | 48%
|
|************************* | 51%
|
|*************************** | 53%
|
|**************************** | 56%
|
|***************************** | 59%
|
|******************************* | 61%
|
|******************************** | 64%
|
|********************************* | 67%
|
|*********************************** | 69%
|
|************************************ | 72%
|
|************************************* | 75%
|
|*************************************** | 77%
|
|**************************************** | 80%
|
|***************************************** | 83%
|
|******************************************* | 85%
|
|******************************************** | 88%
|
|********************************************* | 91%
|
|*********************************************** | 93%
|
|************************************************ | 96%
|
|************************************************* | 99%
|
|**************************************************| 100%
Check the summary and diagnostic plots for the LVM.
summary(lvm.zostera)
$call
boral.default(y = zoo.abnd.flt.zostera, family = "negative.binomial",
lv.control = list(num.lv = 2), row.eff = "fixed", row.ids = matrix(rep(1:5,
times = c(8, 8, 4, 8, 4)), ncol = 1))
$coefficients
coefficients
cols beta0 theta1 theta2 Dispersion
Abra alba 2.198 1.988 0.000 0.805
Abra sp. -1.608 1.860 0.689 20.217
Actiniaria -2.575 1.839 -0.433 16.301
Alitta succinea 1.287 2.962 -0.273 5.909
Ampelisca diadema 3.916 0.203 -2.458 1.892
Amphibalanus improvisus 1.243 2.972 0.827 1.134
Ampithoe sp. 0.444 0.910 -1.998 10.014
Anadara kagoshimensis -2.476 1.031 2.128 15.707
Apherusa bispinosa -2.428 -1.186 0.307 15.831
Apseudopsis ostroumovi 0.521 -5.653 -2.381 2.029
Bittium reticulatum 5.338 -2.188 2.471 0.846
Brachynotus sexdentatus -1.433 0.672 -0.198 14.885
Capitella capitata 1.076 -5.388 -0.928 1.170
Capitella minima 4.479 0.928 -2.812 0.826
Chamelea gallina -0.328 -3.872 0.432 0.531
Chironomidae larvae -2.276 -2.524 -2.237 11.821
Cumella limicola 0.896 -4.777 -0.004 0.436
Cumella pygmaea -1.927 -2.436 1.321 12.821
Cytharella costulata -0.069 -0.433 -0.632 1.365
Diogenes pugilator -0.177 0.968 -2.015 0.423
Eteone flava -2.885 1.639 -2.341 15.594
Eunice vittata -0.779 1.005 -1.601 16.372
Eurydice dollfusi -1.465 -2.009 0.879 15.626
Exogone naidina -0.462 0.723 -3.254 15.923
Gastrosaccus sanctus -2.368 -1.360 0.508 16.826
Genetyllis tuberculata 0.293 0.005 0.259 10.334
Glycera sp. -0.740 2.095 0.024 5.525
Glycera tridactyla -1.041 1.003 2.501 12.689
Glycera unicornis -2.427 0.849 1.844 15.824
Harmothoe imbricata -0.653 0.192 3.196 5.920
Harmothoe reticulata 0.462 0.415 -0.247 0.958
Heteromastus filiformis 3.945 0.050 0.370 0.625
Hirudinea -1.133 0.183 0.720 11.083
Hydrobia acuta -1.418 2.778 -0.307 21.032
Hydrobia sp. -0.577 2.065 -0.426 18.168
Iphinoe tenella 0.352 -0.093 0.673 6.859
Kellia suborbicularis -0.363 -3.112 4.065 3.194
Lagis koreni 1.062 -0.469 3.133 0.278
Leiochone leiopygos 0.318 -1.764 3.210 0.498
Lentidium mediterraneum -2.097 -2.122 0.080 11.619
Lepidochitona cinerea -2.443 -1.418 0.273 15.545
Loripes orbiculatus 1.922 -2.101 0.408 0.397
Lucinella divaricata -1.086 -2.794 0.848 19.669
Magelona papillicornis -2.592 -1.498 -0.331 16.872
Maldane glebifex -2.507 1.761 -0.177 16.326
Melinna palmata 0.652 2.694 -1.236 2.709
Microdeutopus gryllotalpa 1.941 0.402 -2.480 0.546
Micromaldane ornithochaeta -0.860 -0.847 -1.357 15.169
Micronephthys stammeri -1.681 0.374 -1.377 15.521
Microphthalmus fragilis -1.199 -0.311 2.808 21.538
Microphthalmus sp. -1.678 -4.063 -2.984 5.775
Monocorophium acherusicum 2.036 2.349 -3.210 1.058
Mytilaster lineatus 3.458 -2.877 -0.639 1.990
Mytilus galloprovincialis -2.440 1.172 2.084 16.400
Nemertea 1.885 -0.349 1.128 0.641
Nephtys cirrosa -1.858 -1.029 1.939 11.741
Nephtys kersivalensis -2.406 0.799 2.018 15.678
Nereis perivisceralis -2.437 -1.341 0.676 16.179
Nereis pulsatoria -0.673 2.933 1.366 20.309
Nototropis guttatus 0.552 0.331 -1.833 6.902
Oligochaeta 4.784 -1.839 0.190 1.090
Paradoneis harpagonea -2.676 2.076 -0.391 15.422
Parthenina interstincta -0.088 0.859 3.071 13.948
Parvicardium exiguum 2.387 0.840 -0.939 1.586
Perinereis cultrifera 0.507 -0.113 -1.661 6.654
Perioculodes longimanus -0.355 1.073 -2.670 1.471
Phoronida 0.418 -0.105 0.769 0.983
Phyllodoce sp. -3.020 1.597 -2.444 13.608
Platyhelminthes 0.136 2.049 1.017 5.268
Platynereis dumerilii 1.160 0.818 -3.257 0.534
Polititapes aureus -0.311 -0.003 0.208 8.166
Polychaeta larvae -2.476 1.324 -2.128 17.367
Polydora ciliata 3.321 2.466 -0.944 3.736
Polygordius neapolitanus -2.562 1.460 -2.102 16.578
Prionospio cirrifera 3.661 1.111 -0.632 3.972
Protodorvillea kefersteini 1.712 -2.163 0.213 2.182
Pseudocuma longicorne -1.187 0.060 -0.649 12.114
Rissoa membranacea 3.074 0.352 -1.468 1.252
Rissoa splendida -1.177 -1.453 -3.972 2.213
Salvatoria clavata -1.107 -2.064 -4.812 0.653
Schistomeringos rudolphi 0.407 -1.532 4.380 5.729
Sphaerosyllis hystrix -0.252 -3.320 0.508 0.732
Spio filicornis 1.884 2.108 -2.790 3.201
Spisula subtruncata -1.821 -2.397 -0.363 6.609
Stenosoma capito 1.030 -1.190 -1.260 2.092
Syllis gracilis 1.417 -2.385 0.383 2.745
Syllis hyalina -1.882 -2.749 -2.334 19.788
Tellina tenuis -1.540 -1.725 1.252 8.397
Thracia phaseolina -2.348 -1.154 0.295 15.251
Tricolia pullus -1.550 -3.068 -2.264 1.143
Tritia neritea -0.969 1.676 -1.193 14.442
Tritia reticulata 0.147 0.035 1.052 5.493
Turbellaria 0.393 0.927 -0.384 14.925
Upogebia pusilla -2.292 -2.751 -2.773 12.546
$lvs
lv
rows lv1 lv2
1 0.633 -0.251
2 0.719 -0.167
3 0.561 -0.427
4 0.552 -0.125
5 0.310 0.418
6 0.247 0.446
7 0.329 0.500
8 0.410 0.269
9 0.344 -0.052
10 0.306 0.058
11 -0.033 -0.311
12 0.335 0.063
13 -0.032 0.438
14 0.220 0.338
15 0.073 0.375
16 -0.139 0.245
17 0.339 -0.806
18 0.409 -0.803
19 0.595 -1.032
20 0.430 -1.046
21 -0.544 -0.201
22 -0.397 0.003
23 -0.412 -0.046
24 -0.254 0.022
25 -0.341 -0.059
26 -0.433 0.069
27 -0.508 0.140
28 -0.513 0.142
29 -0.474 -0.627
30 -0.555 -0.633
31 -0.342 -0.333
32 -0.474 -0.433
$row.coefficients
$row.coefficients[[1]]
1 2 3 4 5
-1.561 -0.963 -2.739 -0.965 -1.000
$est
[1] "median"
$calc.ics
[1] FALSE
$trial.size
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[44] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[87] 0 0 0 0 0 0 0 0
$num.ord.levels
[1] 0
$prior.control
$prior.control$ssvs.index
[1] -1
attr(,"class")
[1] "summary.boral"
## model fit diagnostic plots
plot(lvm.zostera)
NULL
The residuals plots look fine (no patterns in the residuals vs fitted, so variance is homogeneous, the quantile plot shows a (more or less) normal distribution of the residuals) - the model fits the data pretty well.
Save the zostera LVM.
write_rds(lvm.zostera,
here(save.dir, "lvm_zostera.RDS"))
Examine the biplot obtained by fitting the LVM, as well as the 20 most “important” species.
lvsplot(lvm.zostera, jitter = T, biplot = TRUE, ind.spp = 20)
Only the first 20 ``most important'' latent variable coefficients included in biplot
All in all, the final result resembles the nMDS ordination very much - same stretched clusters (Poda + Otmanli, Vromos pretty much apart, Gradina +- Ropotamo). I don’t see much difference with the nMDS. The main difference seems to be the distance between the 2 years for Poda ana Otmanli - the LVM enlarges it. Have to remember to test for year effect! The run time is actually not that bad for the seagrasses. The species singled out as significant are probably somewhat different - have to check!
Redo the biplot, because this one is not very pretty. I’m not adding the species on top, first because I’m too lazy to figure out the procedure for ordering them, and second because the plot gets too busy.
## extract the LV coordinates of the stations from the model, so that the plot can be redone in ggplot
lvs.coord.zostera <- as_tibble(lvm.zostera$lv.median)
## add the stations from the original zoobenthic table (order was not modified)
(lvs.coord.zostera <- lvs.coord.zostera %>%
bind_cols(zoo.abnd.zostera %>% select(station))
)
Make the plot and save it.
(plot.lvm.zostera <- ggplot(lvs.coord.zostera) +
geom_point(aes(x = lv1, y = lv2, colour = station)) +
scale_color_brewer(palette = "Set2", name = "station",
labels = paste0("Z", as.numeric((unique(lvs.coord.zostera %>% pull(station)))))) +
labs(x = "LV1", y = "LV2")
)
ggsave(file = here(figures.dir, "lvm_zostera.png"),
plot.lvm.zostera,
width = 15, units = "cm", dpi = 300)
Saving 15 x 11.4 cm image
Well, this is a weird one - this plot is flipped around 0 compared to the one that boral’s plotting function gives. Otherwise nothing changes - the spatial relationships between samples are preserved. I suppose it doesn’t matter much - the axes are arbitrary after all, but strange that it happens.
## save the LVM plot for the seagrass
ggsave(file = here(figures.dir, "lvm_zostera.png"),
plot.lvm.zostera,
width = 15, units = "cm", dpi = 300)
Saving 15 x 17.8 cm image
Fit GLMs to the sites x species matrix to try and explain the observed differences in community structure by the variation of the environmental parameters.
These functions all come from package mvabund.
Import the environmental data - the one cleaned, prepared and saved in the previous notebook (classical multivariate methods). It contains long-term averages for the water column data (as long-term as available, at least) at each station, repeated for each replicate, the sediment data (2013-2014), and the seagrass data (2013-2014), again repeated to the same number of replicates. Only the variables determined to be significant by PCA are kept.
env.zostera <- read_csv(here(save.dir, "env_data_ordinations_zostera.csv"))
Parsed with column specification:
cols(
station = [31mcol_character()[39m,
year = [32mcol_double()[39m,
Ntotal = [32mcol_double()[39m,
chl_a = [32mcol_double()[39m,
secchi = [32mcol_double()[39m,
LUSI = [32mcol_double()[39m,
TOM = [32mcol_double()[39m,
moisture_content = [32mcol_double()[39m,
mean_grain_size = [32mcol_double()[39m,
sand = [32mcol_double()[39m,
silt_clay = [32mcol_double()[39m,
shoot_count = [32mcol_double()[39m,
ag_biomass_wet = [32mcol_double()[39m,
bg_biomass_wet = [32mcol_double()[39m
)
## convert station to factor
(env.zostera <- env.zostera %>%
mutate(station = factor(station,
levels = c("Poda", "Otmanli", "Vromos", "Gradina", "Ropotamo")))
)
Station is a factor, the rest of the variables are numeric.
Turn the zoobenthic data (minus the all-0 taxa) into a matrix - easier for the mvabund package and methods to deal with.
## there is already one subset of filtered count data (32 x 94) - use it
zoo.mvabnd.zostera <- mvabund(zoo.abnd.flt.zostera)
First, let’s see if the groups from the latent variable model (more or less equal to the clusters from the classical ordination) are valid, and which species exhibit a response.
I’m going to try something new here - 1) loose clusters from the LVM ordination, 1 = Poda-Otmanli, 2 = Vromos, 3 = Gradina-Ropotamo. 2) stations as clusters, as I did before for the seagrass data, although I don’t believe it’s valid/justified to do so… 3) another configuration of clusters from the LVM ordination: 1 = Z1-Z2, 2 = Z3, 3 = Z4, 4 = Z5.
## construct the vectors of the clusters by hand - first, situation 1 above
lvm.clusters.zostera.1 <- rep(1:3, times = c(16, 4, 12))
(lvm.clusters.zostera.1 <- factor(lvm.clusters.zostera.1))
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3
Levels: 1 2 3
## again, for case 2
lvm.clusters.zostera.2 <- rep(1:5, times = c(8, 8, 4, 8, 4))
(lvm.clusters.zostera.2 <- factor(lvm.clusters.zostera.2))
[1] 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4 5 5 5 5
Levels: 1 2 3 4 5
## again, for case 3
lvm.clusters.zostera.3 <- rep(1:4, times = c(16, 4, 8, 4))
(lvm.clusters.zostera.3 <- factor(lvm.clusters.zostera.3))
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4
Levels: 1 2 3 4
LVM clusters - case 1 Check the model assumptions. 1. Mean-variance assumption => determines the choice of family parameter. Can be checked by plotting residuals vs fits: if little pattern - the chosen mean-variance assumption is plausible.
Another way: direct plotting (variance ~ mean), for each species within each factor level.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1, table = TRUE)
START SECTION 2
Plotting if overlay is TRUE
using grouping variable lvm.clusters.zostera.1 105 mean values were 0 and could
not be included in the log-plot
using grouping variable lvm.clusters.zostera.1 105 variance values were 0 and could not
be included in the log-plot
FINISHED SECTION 2
$mean
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 57.6875 50.8125 30.0625 6.25000 3.25000
2 0.0000 88.0000 0.7500 80.50000 2.25000
3 183.6667 56.7500 99.5000 23.08333 51.83333
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata Monocorophium.acherusicum
1 17.81250 25.375000 26.562500 3.937500
2 2.25000 0.000000 0.000000 78.250000
3 17.08333 6.666667 3.916667 1.416667
Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi Spio.filicornis
1 7.437500 0.43750 0.06250 1.4375000
2 9.500000 0.00000 0.00000 38.2500000
3 8.083333 19.83333 18.66667 0.8333333
Microdeutopus.gryllotalpa Abra.alba Cumella.limicola Loripes.orbiculatus
1 2.625000 6.312500 0.375 1.125000
2 4.250000 0.000000 0.000 0.500000
3 5.916667 1.083333 9.000 7.333333
Parvicardium.exiguum Protodorvillea.kefersteini Platynereis.dumerilii Nemertea
1 2.812500 0.75 1.562500 2.375000
2 5.500000 0.00 9.250000 0.500000
3 2.416667 7.00 2.333333 2.583333
Syllis.gracilis Alitta.succinea Amphibalanus.improvisus Stenosoma.capito Lagis.koreni
1 1.625 3.5625000 3.3125000 1.500000 2.1250000
2 0.000 0.0000000 0.2500000 0.000000 0.0000000
3 3.000 0.3333333 0.3333333 1.833333 0.9166667
Schistomeringos.rudolphi Salvatoria.clavata Leiochone.leiopygos Melinna.palmata
1 1.6875000 0.000000 0.625 0.875
2 0.0000000 1.500000 0.000 2.750
3 0.8333333 2.333333 1.250 0.000
Microphthalmus.sp. Kellia.suborbicularis Nototropis.guttatus Chamelea.gallina
1 0.000000 0.312500 0.7500000 0.0625
2 0.000000 0.000000 0.0000000 0.0000
3 2.083333 1.583333 0.9166667 1.7500
Perinereis.cultrifera Sphaerosyllis.hystrix Ampithoe.sp. Harmothoe.reticulata Phoronida
1 0.437500 0.187500 0.1875000 0.7500000 0.5625000
2 0.000000 0.000000 2.7500000 0.0000000 0.0000000
3 1.083333 1.333333 0.3333333 0.4166667 0.5833333
Perioculodes.longimanus Rissoa.splendida Diogenes.pugilator Iphinoe.tenella
1 0.3750000 0.0000000 0.3750000 0.6875000
2 1.2500000 1.0000000 0.7500000 0.0000000
3 0.3333333 0.9166667 0.4166667 0.1666667
Platyhelminthes Exogone.naidina Genetyllis.tuberculata Parthenina.interstincta
1 0.6875000 0.0625000 0.50 0.6250000
2 0.0000000 2.2500000 0.25 0.0000000
3 0.1666667 0.1666667 0.25 0.1666667
Tritia.reticulata Cytharella.costulata Syllis.hyalina Tricolia.pullus Turbellaria
1 0.5625 0.1875000 0.0000000 0.0625 0.375
2 0.0000 0.2500000 0.0000000 0.0000 0.250
3 0.2500 0.5833333 0.8333333 0.7500 0.250
Harmothoe.imbricata Hydrobia.sp. Lucinella.divaricata Nereis.pulsatoria
1 0.3125000 0.37500000 0.0000000 0.4375
2 0.0000000 0.00000000 0.0000000 0.0000
3 0.1666667 0.08333333 0.5833333 0.0000
Polititapes.aureus Upogebia.pusilla Glycera.sp. Microphthalmus.fragilis Eunice.vittata
1 0.37500000 0.0000000 0.375 0.375 0.1875000
2 0.00000000 0.0000000 0.000 0.000 0.0000000
3 0.08333333 0.5833333 0.000 0.000 0.1666667
Glycera.tridactyla Tritia.neritea Chironomidae.larvae Hydrobia.acuta
1 0.3125 0.25000000 0.0000000 0.25
2 0.0000 0.00000000 0.0000000 0.00
3 0.0000 0.08333333 0.3333333 0.00
Micromaldane.ornithochaeta Cumella.pygmaea Eurydice.dollfusi Hirudinea
1 0.00 0.00 0.00 0.12500000
2 0.25 0.00 0.00 0.00000000
3 0.25 0.25 0.25 0.08333333
Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis Brachynotus.sexdentatus
1 0.0625000 0.00 0.0625000 0.06250000
2 0.0000000 0.00 0.0000000 0.00000000
3 0.1666667 0.25 0.1666667 0.08333333
Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa Abra.sp. Actiniaria
1 0.0000000 0.06250000 0.06250000 0.00 0.0625
2 0.0000000 0.00000000 0.00000000 0.25 0.0000
3 0.1666667 0.08333333 0.08333333 0.00 0.0000
Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava Gastrosaccus.sanctus
1 0.0625 0.00000000 0.00 0.00000000
2 0.0000 0.00000000 0.25 0.00000000
3 0.0000 0.08333333 0.00 0.08333333
Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis Maldane.glebifex
1 0.0625 0.00000000 0.00000000 0.0625
2 0.0000 0.00000000 0.00000000 0.0000
3 0.0000 0.08333333 0.08333333 0.0000
Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.0625 0.0625 0.00000000
2 0.0000 0.0000 0.00000000
3 0.0000 0.0000 0.08333333
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.0625 0.00 0.00 0.00
2 0.0000 0.25 0.25 0.25
3 0.0000 0.00 0.00 0.00
Thracia.phaseolina
1 0.00000000
2 0.00000000
3 0.08333333
$var
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 2642.229 2899.229 700.8625 45.53333 35.000000
2 0.000 3094.000 2.2500 395.00000 1.583333
3 23448.242 3601.295 13116.2727 497.53788 1873.060606
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata Monocorophium.acherusicum
1 177.495833 919.5833 1000.26250 17.92917
2 1.583333 0.0000 0.00000 3980.91667
3 171.901515 146.9697 38.81061 2.44697
Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi Spio.filicornis
1 271.72917 1.595833 0.0625 2.529167
2 11.66667 0.000000 0.0000 306.250000
3 71.71970 423.606061 537.1515 2.515152
Microdeutopus.gryllotalpa Abra.alba Cumella.limicola Loripes.orbiculatus
1 10.78333 24.362500 1.183333 2.383333
2 16.25000 0.000000 0.000000 1.000000
3 48.44697 1.356061 29.090909 16.787879
Parvicardium.exiguum Protodorvillea.kefersteini Platynereis.dumerilii Nemertea
1 18.829167 2.60000 6.929167 12.2500000
2 3.666667 0.00000 32.916667 0.3333333
3 6.265152 29.27273 5.696970 4.0833333
Syllis.gracilis Alitta.succinea Amphibalanus.improvisus Stenosoma.capito Lagis.koreni
1 7.45 30.1291667 23.4291667 19.733333 6.5166667
2 0.00 0.0000000 0.2500000 0.000000 0.0000000
3 18.00 0.4242424 0.4242424 6.151515 0.8106061
Schistomeringos.rudolphi Salvatoria.clavata Leiochone.leiopygos Melinna.palmata
1 11.29583 0.0000000 0.650000 1.183333
2 0.00000 0.3333333 0.000000 0.250000
3 2.69697 13.1515152 3.295455 0.000000
Microphthalmus.sp. Kellia.suborbicularis Nototropis.guttatus Chamelea.gallina
1 0.00000 0.6291667 1.666667 0.062500
2 0.00000 0.0000000 0.000000 0.000000
3 12.62879 10.0833333 5.356061 1.659091
Perinereis.cultrifera Sphaerosyllis.hystrix Ampithoe.sp. Harmothoe.reticulata Phoronida
1 2.262500 0.2958333 0.2958333 0.7333333 0.5291667
2 0.000000 0.0000000 11.5833333 0.0000000 0.0000000
3 1.901515 1.6969697 0.4242424 0.8106061 0.8106061
Perioculodes.longimanus Rissoa.splendida Diogenes.pugilator Iphinoe.tenella
1 0.5166667 0.000000 0.3833333 2.2291667
2 3.5833333 2.000000 0.2500000 0.0000000
3 0.4242424 2.992424 0.4469697 0.3333333
Platyhelminthes Exogone.naidina Genetyllis.tuberculata Parthenina.interstincta
1 1.4291667 0.0625000 2.2666667 2.6500000
2 0.0000000 8.2500000 0.2500000 0.0000000
3 0.1515152 0.3333333 0.3863636 0.3333333
Tritia.reticulata Cytharella.costulata Syllis.hyalina Tricolia.pullus Turbellaria
1 1.4625000 0.1625000 0.000000 0.062500 1.05
2 0.0000000 0.2500000 0.000000 0.000000 0.25
3 0.2045455 0.6287879 8.333333 1.477273 0.75
Harmothoe.imbricata Hydrobia.sp. Lucinella.divaricata Nereis.pulsatoria
1 0.6291667 1.18333333 0.000000 1.4625
2 0.0000000 0.00000000 0.000000 0.0000
3 0.1515152 0.08333333 2.265152 0.0000
Polititapes.aureus Upogebia.pusilla Glycera.sp. Microphthalmus.fragilis Eunice.vittata
1 0.65000000 0.000000 0.3833333 2.25 0.2958333
2 0.00000000 0.000000 0.0000000 0.00 0.0000000
3 0.08333333 2.265152 0.0000000 0.00 0.3333333
Glycera.tridactyla Tritia.neritea Chironomidae.larvae Hydrobia.acuta
1 0.6291667 0.60000000 0.0000000 1
2 0.0000000 0.00000000 0.0000000 0
3 0.0000000 0.08333333 0.7878788 0
Micromaldane.ornithochaeta Cumella.pygmaea Eurydice.dollfusi Hirudinea
1 0.0000000 0.0000000 0.0000000 0.11666667
2 0.2500000 0.0000000 0.0000000 0.00000000
3 0.3863636 0.3863636 0.3863636 0.08333333
Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis Brachynotus.sexdentatus
1 0.0625000 0.0000000 0.0625000 0.06250000
2 0.0000000 0.0000000 0.0000000 0.00000000
3 0.1515152 0.2045455 0.1515152 0.08333333
Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa Abra.sp. Actiniaria
1 0.0000000 0.06250000 0.06250000 0.00 0.0625
2 0.0000000 0.00000000 0.00000000 0.25 0.0000
3 0.1515152 0.08333333 0.08333333 0.00 0.0000
Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava Gastrosaccus.sanctus
1 0.0625 0.00000000 0.00 0.00000000
2 0.0000 0.00000000 0.25 0.00000000
3 0.0000 0.08333333 0.00 0.08333333
Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis Maldane.glebifex
1 0.0625 0.00000000 0.00000000 0.0625
2 0.0000 0.00000000 0.00000000 0.0000
3 0.0000 0.08333333 0.08333333 0.0000
Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.0625 0.0625 0.00000000
2 0.0000 0.0000 0.00000000
3 0.0000 0.0000 0.08333333
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.0625 0.00 0.00 0.00
2 0.0000 0.25 0.25 0.25
3 0.0000 0.00 0.00 0.00
Thracia.phaseolina
1 0.00000000
2 0.00000000
3 0.08333333
It’s not perfect, but it’s not too terrible either.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.1 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.1,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.1)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.1, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(1:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(1, 1))
I really don’t like the rainbow palette, but I would like to include these plots in my thesis results.. Will have to do something about it, just not right now.
Save the model!
write_rds(glms.lvm.zostera.1,
here(save.dir, "glms_lvm_zostera_1.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.1.summary <- summary(glms.lvm.zostera.1,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
Resampling run 0 finished. Time elapsed: 0.00 min ...
Resampling run 100 finished. Time elapsed: 0.28 min ...
Resampling run 200 finished. Time elapsed: 0.57 min ...
Resampling run 300 finished. Time elapsed: 0.85 min ...
Resampling run 400 finished. Time elapsed: 1.14 min ...
Resampling run 500 finished. Time elapsed: 1.46 min ...
Resampling run 600 finished. Time elapsed: 1.74 min ...
Resampling run 700 finished. Time elapsed: 2.05 min ...
Resampling run 800 finished. Time elapsed: 2.39 min ...
Resampling run 900 finished. Time elapsed: 2.70 min ...
Time elapsed: 0 hr 2 min 58 sec
Test statistics:
LR value Pr(>LR)
(Intercept) 1263.9 0.001 ***
lvm.clusters.zostera.12 319.1 0.001 ***
lvm.clusters.zostera.13 386.9 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate test statistic:
(Intercept) lvm.clusters.zostera.12
LR value Pr(>LR) LR value Pr(>LR)
Abra.alba 46.653 0.001 17.680 0.002
Abra.sp. 6.107 0.514 3.219 0.648
Actiniaria 2.833 0.952 0.446 0.989
Alitta.succinea 7.171 0.397 4.935 0.336
Ampelisca.diadema 45.210 0.001 20.026 0.001
Amphibalanus.improvisus 15.311 0.049 5.964 0.286
Ampithoe.sp. 5.196 0.632 5.826 0.293
Anadara.kagoshimensis 2.833 0.952 0.446 0.989
Apherusa.bispinosa 4.442 0.795 0.000 0.998
Apseudopsis.ostroumovi 11.749 0.137 0.426 0.989
Bittium.reticulatum 156.907 0.001 24.100 0.001
Brachynotus.sexdentatus 4.472 0.795 0.446 0.989
Capitella.capitata 2.801 0.952 2.426 0.768
Capitella.minima 153.136 0.001 0.735 0.976
Chamelea.gallina 20.892 0.020 0.446 0.989
Chironomidae.larvae 5.475 0.592 0.000 1.000
Cumella.limicola 6.114 0.514 2.513 0.755
Cumella.pygmaea 6.270 0.503 0.000 0.998
Cytharella.costulata 10.234 0.152 0.059 0.998
Diogenes.pugilator 7.396 0.378 0.877 0.974
Eteone.flava 6.107 0.514 3.219 0.648
Eunice.vittata 2.143 0.987 0.865 0.974
Eurydice.dollfusi 6.270 0.503 0.000 0.998
Exogone.naidina 5.740 0.564 5.438 0.297
Gastrosaccus.sanctus 4.442 0.795 0.000 0.998
Genetyllis.tuberculata 0.906 0.996 0.171 0.998
Glycera.sp. 4.780 0.737 2.664 0.736
Glycera.tridactyla 2.091 0.987 1.471 0.921
Glycera.unicornis 2.833 0.952 0.446 0.989
Harmothoe.imbricata 2.901 0.947 1.687 0.879
Harmothoe.reticulata 0.846 0.996 4.747 0.360
Heteromastus.filiformis 117.049 0.001 10.036 0.049
Hirudinea 5.554 0.582 0.893 0.974
Hydrobia.acuta 0.471 0.998 0.467 0.989
Hydrobia.sp. 0.891 0.996 1.130 0.951
Iphinoe.tenella 0.518 0.998 2.798 0.705
Kellia.suborbicularis 2.134 0.987 1.390 0.921
Lagis.koreni 6.465 0.495 8.431 0.097
Leiochone.leiopygos 1.439 0.992 3.585 0.559
Lentidium.mediterraneum 8.262 0.263 0.000 1.000
Lepidochitona.cinerea 4.442 0.795 0.000 0.998
Loripes.orbiculatus 0.187 0.998 1.185 0.939
Lucinella.divaricata 4.632 0.770 0.000 0.998
Magelona.papillicornis 4.442 0.795 0.000 0.998
Maldane.glebifex 2.833 0.952 0.446 0.989
Melinna.palmata 0.261 0.998 6.839 0.191
Microdeutopus.gryllotalpa 9.960 0.163 0.553 0.989
Micromaldane.ornithochaeta 8.472 0.246 2.898 0.672
Micronephthys.stammeri 4.472 0.795 0.446 0.989
Microphthalmus.fragilis 0.234 0.998 0.467 0.989
Microphthalmus.sp. 7.752 0.352 0.000 0.998
Monocorophium.acherusicum 26.311 0.010 26.112 0.001
Mytilaster.lineatus 12.543 0.116 0.198 0.998
Mytilus.galloprovincialis 2.833 0.952 0.446 0.989
Nemertea 8.792 0.238 3.214 0.648
Nephtys.cirrosa 4.472 0.795 0.446 0.989
Nephtys.kersivalensis 2.833 0.952 0.446 0.989
Nereis.perivisceralis 4.442 0.795 0.000 0.998
Nereis.pulsatoria 0.491 0.998 0.980 0.968
Nototropis.guttatus 0.217 0.998 2.401 0.774
Oligochaeta 121.705 0.001 13.320 0.009
Paradoneis.harpagonea 2.833 0.952 0.446 0.989
Parthenina.interstincta 0.277 0.998 1.459 0.921
Parvicardium.exiguum 13.759 0.086 1.460 0.921
Perinereis.cultrifera 1.679 0.992 1.968 0.853
Perioculodes.longimanus 3.532 0.892 1.819 0.879
Phoronida 2.975 0.947 3.857 0.496
Phyllodoce.sp. 6.107 0.514 3.219 0.648
Platyhelminthes 0.614 0.998 3.040 0.648
Platynereis.dumerilii 1.732 0.992 7.023 0.183
Polititapes.aureus 2.721 0.954 2.061 0.839
Polychaeta.larvae 6.107 0.514 3.219 0.648
Polydora.ciliata 69.320 0.001 11.432 0.027
Polygordius.neapolitanus 6.107 0.514 3.219 0.648
Prionospio.cirrifera 67.239 0.001 11.104 0.029
Protodorvillea.kefersteini 0.496 0.998 3.781 0.508
Pseudocuma.longicorne 6.930 0.445 0.446 0.989
Rissoa.membranacea 49.479 0.001 0.153 0.998
Rissoa.splendida 13.625 0.088 8.728 0.086
Salvatoria.clavata 14.949 0.056 11.446 0.027
Schistomeringos.rudolphi 0.600 0.998 2.474 0.761
Sphaerosyllis.hystrix 10.204 0.153 1.307 0.921
Spio.filicornis 1.529 0.992 21.365 0.001
Spisula.subtruncata 11.731 0.137 0.000 1.000
Stenosoma.capito 0.912 0.996 4.619 0.371
Syllis.gracilis 1.237 0.996 4.542 0.389
Syllis.hyalina 1.975 0.987 0.000 1.000
Tellina.tenuis 6.930 0.445 0.446 0.989
Thracia.phaseolina 4.442 0.795 0.000 0.998
Tricolia.pullus 9.618 0.176 0.432 0.989
Tritia.neritea 1.843 0.991 1.073 0.955
Tritia.reticulata 1.128 0.996 2.531 0.754
Turbellaria 1.295 0.996 0.051 0.998
Upogebia.pusilla 4.632 0.770 0.000 1.000
lvm.clusters.zostera.13
LR value Pr(>LR)
Abra.alba 15.416 0.010
Abra.sp. 0.000 1.000
Actiniaria 1.119 1.000
Alitta.succinea 5.293 0.616
Ampelisca.diadema 9.548 0.133
Amphibalanus.improvisus 12.041 0.053
Ampithoe.sp. 0.374 1.000
Anadara.kagoshimensis 1.119 1.000
Apherusa.bispinosa 1.695 1.000
Apseudopsis.ostroumovi 26.735 0.001
Bittium.reticulatum 7.624 0.255
Brachynotus.sexdentatus 0.041 1.000
Capitella.capitata 21.010 0.003
Capitella.minima 0.058 1.000
Chamelea.gallina 22.405 0.001
Chironomidae.larvae 3.638 0.907
Cumella.limicola 26.901 0.001
Cumella.pygmaea 3.669 0.904
Cytharella.costulata 2.997 0.965
Diogenes.pugilator 0.030 1.000
Eteone.flava 0.000 1.000
Eunice.vittata 0.008 1.000
Eurydice.dollfusi 3.669 0.904
Exogone.naidina 0.518 1.000
Gastrosaccus.sanctus 1.695 1.000
Genetyllis.tuberculata 0.391 1.000
Glycera.sp. 6.512 0.407
Glycera.tridactyla 3.637 0.907
Glycera.unicornis 1.119 1.000
Harmothoe.imbricata 0.391 1.000
Harmothoe.reticulata 1.102 1.000
Heteromastus.filiformis 0.017 1.000
Hirudinea 0.114 1.000
Hydrobia.acuta 1.163 1.000
Hydrobia.sp. 0.826 1.000
Iphinoe.tenella 2.093 0.991
Kellia.suborbicularis 2.456 0.983
Lagis.koreni 2.825 0.967
Leiochone.leiopygos 1.600 1.000
Lentidium.mediterraneum 3.389 0.934
Lepidochitona.cinerea 1.695 1.000
Loripes.orbiculatus 21.059 0.003
Lucinella.divaricata 3.696 0.902
Magelona.papillicornis 1.695 1.000
Maldane.glebifex 1.119 1.000
Melinna.palmata 15.669 0.009
Microdeutopus.gryllotalpa 3.082 0.961
Micromaldane.ornithochaeta 4.262 0.794
Micronephthys.stammeri 0.041 1.000
Microphthalmus.fragilis 1.163 1.000
Microphthalmus.sp. 8.208 0.201
Monocorophium.acherusicum 5.392 0.602
Mytilaster.lineatus 18.655 0.003
Mytilus.galloprovincialis 1.119 1.000
Nemertea 0.038 1.000
Nephtys.cirrosa 0.041 1.000
Nephtys.kersivalensis 1.119 1.000
Nereis.perivisceralis 1.695 1.000
Nereis.pulsatoria 2.423 0.984
Nototropis.guttatus 0.050 1.000
Oligochaeta 7.645 0.252
Paradoneis.harpagonea 1.119 1.000
Parthenina.interstincta 0.797 1.000
Parvicardium.exiguum 0.131 1.000
Perinereis.cultrifera 1.135 1.000
Perioculodes.longimanus 0.023 1.000
Phoronida 0.005 1.000
Phyllodoce.sp. 0.000 1.000
Platyhelminthes 2.344 0.987
Platynereis.dumerilii 0.630 1.000
Polititapes.aureus 1.865 0.998
Polychaeta.larvae 0.000 1.000
Polydora.ciliata 6.959 0.357
Polygordius.neapolitanus 0.000 1.000
Prionospio.cirrifera 3.829 0.882
Protodorvillea.kefersteini 11.945 0.055
Pseudocuma.longicorne 0.689 1.000
Rissoa.membranacea 0.036 1.000
Rissoa.splendida 10.951 0.077
Salvatoria.clavata 16.907 0.003
Schistomeringos.rudolphi 0.375 1.000
Sphaerosyllis.hystrix 9.295 0.147
Spio.filicornis 1.261 1.000
Spisula.subtruncata 5.084 0.644
Stenosoma.capito 0.093 1.000
Syllis.gracilis 0.821 1.000
Syllis.hyalina 1.773 1.000
Tellina.tenuis 0.689 1.000
Thracia.phaseolina 1.695 1.000
Tricolia.pullus 6.175 0.457
Tritia.neritea 0.559 1.000
Tritia.reticulata 0.825 1.000
Turbellaria 0.114 1.000
Upogebia.pusilla 3.696 0.902
Arguments: with 999 resampling iterations using pit.trap resampling and response assumed to be uncorrelated
Likelihood Ratio statistic: 703.6, p-value: 0.001
Univariate test statistic:
Abra.alba Abra.sp. Actiniaria Alitta.succinea Ampelisca.diadema
LR value 25.931 4.159 1.386 8.646 20.059
Pr(>LR) 0.001 0.900 0.998 0.323 0.002
Amphibalanus.improvisus Ampithoe.sp. Anadara.kagoshimensis Apherusa.bispinosa
LR value 14.882 6.321 1.386 1.962
Pr(>LR) 0.030 0.660 0.998 0.994
Apseudopsis.ostroumovi Bittium.reticulatum Brachynotus.sexdentatus
LR value 30.855 31.183 0.575
Pr(>LR) 0.001 0.001 1.000
Capitella.capitata Capitella.minima Chamelea.gallina Chironomidae.larvae
LR value 26.213 0.739 26.230 4.195
Pr(>LR) 0.001 1.000 0.001 0.891
Cumella.limicola Cumella.pygmaea Cytharella.costulata Diogenes.pugilator
LR value 32.545 4.228 3.122 0.896
Pr(>LR) 0.001 0.886 0.969 1.000
Eteone.flava Eunice.vittata Eurydice.dollfusi Exogone.naidina
LR value 4.159 0.871 4.228 5.944
Pr(>LR) 0.900 1.000 0.886 0.694
Gastrosaccus.sanctus Genetyllis.tuberculata Glycera.sp. Glycera.tridactyla
LR value 1.962 0.467 7.992 4.484
Pr(>LR) 0.994 1.000 0.397 0.851
Glycera.unicornis Harmothoe.imbricata Harmothoe.reticulata
LR value 1.386 1.841 5.130
Pr(>LR) 0.998 0.994 0.777
Heteromastus.filiformis Hirudinea Hydrobia.acuta Hydrobia.sp. Iphinoe.tenella
LR value 10.222 0.915 1.438 1.681 4.160
Pr(>LR) 0.194 1.000 0.998 0.998 0.895
Kellia.suborbicularis Lagis.koreni Leiochone.leiopygos Lentidium.mediterraneum
LR value 4.584 9.705 6.218 3.923
Pr(>LR) 0.842 0.226 0.660 0.918
Lepidochitona.cinerea Loripes.orbiculatus Lucinella.divaricata
LR value 1.962 25.259 4.257
Pr(>LR) 0.994 0.001 0.882
Magelona.papillicornis Maldane.glebifex Melinna.palmata
LR value 1.962 1.386 24.203
Pr(>LR) 0.994 0.998 0.001
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta Micronephthys.stammeri
LR value 3.088 4.544 0.575
Pr(>LR) 0.969 0.848 1.000
Microphthalmus.fragilis Microphthalmus.sp. Monocorophium.acherusicum
LR value 1.437 9.396 38.329
Pr(>LR) 0.998 0.260 0.001
Mytilaster.lineatus Mytilus.galloprovincialis Nemertea Nephtys.cirrosa
LR value 21.933 1.386 3.554 0.575
Pr(>LR) 0.001 0.998 0.941 1.000
Nephtys.kersivalensis Nereis.perivisceralis Nereis.pulsatoria
LR value 1.386 1.962 2.987
Pr(>LR) 0.998 0.994 0.969
Nototropis.guttatus Oligochaeta Paradoneis.harpagonea Parthenina.interstincta
LR value 2.610 21.776 1.386 1.965
Pr(>LR) 0.975 0.001 0.998 0.994
Parvicardium.exiguum Perinereis.cultrifera Perioculodes.longimanus Phoronida
LR value 2.086 3.680 2.254 4.088
Pr(>LR) 0.994 0.941 0.994 0.905
Phyllodoce.sp. Platyhelminthes Platynereis.dumerilii Polititapes.aureus
LR value 4.159 4.534 7.280 3.315
Pr(>LR) 0.900 0.848 0.498 0.955
Polychaeta.larvae Polydora.ciliata Polygordius.neapolitanus Prionospio.cirrifera
LR value 4.159 15.690 4.159 13.193
Pr(>LR) 0.900 0.021 0.900 0.058
Protodorvillea.kefersteini Pseudocuma.longicorne Rissoa.membranacea
LR value 17.536 1.491 0.158
Pr(>LR) 0.008 0.998 1.000
Rissoa.splendida Salvatoria.clavata Schistomeringos.rudolphi
LR value 11.549 17.003 2.634
Pr(>LR) 0.111 0.008 0.975
Sphaerosyllis.hystrix Spio.filicornis Spisula.subtruncata Stenosoma.capito
LR value 12.713 28.109 5.885 5.011
Pr(>LR) 0.075 0.001 0.694 0.796
Syllis.gracilis Syllis.hyalina Tellina.tenuis Thracia.phaseolina Tricolia.pullus
LR value 5.956 2.048 1.491 1.962 7.956
Pr(>LR) 0.694 0.994 0.998 0.994 0.401
Tritia.neritea Tritia.reticulata Turbellaria Upogebia.pusilla
LR value 1.406 2.932 0.135 4.257
Pr(>LR) 0.998 0.969 1.000 0.882
Arguments:
Test statistics calculated assuming response assumed to be uncorrelated
P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).
The factor is highly significant according to the models.
This also allows us to see which species exhibit a response to the chosen factor. The LR (likelihood ratio) statistic is used as a measure of the strength of individual taxon contributions to the observed patterns. I’ll save the summary for safekeeping, but I’ll also run an anova - to get an analysis of deviance table on the model fit (also better for extracting the species contributions, or at least I know how to do it).
write_rds(glms.lvm.zostera.1.summary,
here(save.dir, "glms_lvm_zostera_1_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.1.aov <- anova.manyglm(glms.lvm.zostera.1,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.24 minutes...
Resampling run 300 finished. Time elapsed: 0.35 minutes...
Resampling run 400 finished. Time elapsed: 0.47 minutes...
Resampling run 500 finished. Time elapsed: 0.58 minutes...
Resampling run 600 finished. Time elapsed: 0.70 minutes...
Resampling run 700 finished. Time elapsed: 0.82 minutes...
Resampling run 800 finished. Time elapsed: 0.93 minutes...
Resampling run 900 finished. Time elapsed: 1.05 minutes...
Time elapsed: 0 hr 1 min 9 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.1,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.1 29 2 703.6 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 25.931 0.001 4.159 0.983 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 8.646 0.433 20.059 0.006
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 14.882 0.036 6.321 0.799
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 1.962 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 30.855 0.001 31.183 0.001
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.575 1.000 26.213 0.001
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.739 1.000 26.23 0.001
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.195 0.981 32.545 0.001
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.228 0.981 3.122 0.996
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.896 1.000 4.159 0.986 0.871
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 4.228 0.981 5.944 0.838
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.962 1.000 0.467 1.000
Glycera.sp. Glycera.tridactyla Glycera.unicornis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 7.992 0.560 4.484 0.970 1.386
Harmothoe.imbricata Harmothoe.reticulata
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.000 1.841 1.000 5.13
Heteromastus.filiformis Hirudinea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.923 10.222 0.254 0.915 1.000
Hydrobia.acuta Hydrobia.sp. Iphinoe.tenella
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.438 1.000 1.681 1.000 4.16
Kellia.suborbicularis Lagis.koreni
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.981 4.584 0.953 9.705 0.311
Leiochone.leiopygos Lentidium.mediterraneum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 6.218 0.803 3.923 0.986
Lepidochitona.cinerea Loripes.orbiculatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.962 1.000 25.259 0.001
Lucinella.divaricata Magelona.papillicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.257 0.979 1.962 1.000
Maldane.glebifex Melinna.palmata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 24.203 0.001
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 3.088 0.998 4.544
Micronephthys.stammeri Microphthalmus.fragilis
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.965 0.575 1.000 1.437
Microphthalmus.sp. Monocorophium.acherusicum
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.000 9.396 0.348 38.329
Mytilaster.lineatus Mytilus.galloprovincialis
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.001 21.933 0.003 1.386
Nemertea Nephtys.cirrosa
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 3.554 0.992 0.575 1.000
Nephtys.kersivalensis Nereis.perivisceralis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.386 1.000 1.962 1.000
Nereis.pulsatoria Nototropis.guttatus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 2.987 0.999 2.61 1.000
Oligochaeta Paradoneis.harpagonea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 21.776 0.003 1.386 1.000
Parthenina.interstincta Parvicardium.exiguum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.965 1.000 2.086 1.000
Perinereis.cultrifera Perioculodes.longimanus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 3.68 0.992 2.254 1.000
Phoronida Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 4.088 0.986 4.159 0.986 4.534
Platynereis.dumerilii Polititapes.aureus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 0.966 7.28 0.630 3.315
Polychaeta.larvae Polydora.ciliata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.996 4.159 0.982 15.69 0.031
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 4.159 0.982 13.193 0.085
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 17.536 0.010 1.491 1.000
Rissoa.membranacea Rissoa.splendida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 0.158 1.000 11.549 0.153
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 17.003 0.014 2.634 1.000
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 12.713 0.097 28.109 0.001
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 5.885 0.844 5.011 0.923
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 5.956 0.838 2.048 1.000 1.491
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 1.962 1.000 7.956 0.560
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.1 1.406 1.000 2.932 0.999 0.135
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.1 1.000 4.257 0.979
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Save the ANOVA, too.
write_rds(glms.lvm.zostera.1.aov,
here(save.dir, "glms_lvm_zostera_1_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.1 <- top_n_sp_glm(glms.lvm.zostera.1.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.76"
Monocorophium.acherusicum Cumella.limicola Bittium.reticulatum
38.328670 32.545296 31.182981
Apseudopsis.ostroumovi Spio.filicornis Chamelea.gallina
30.854772 28.108956 26.230014
Capitella.capitata Abra.alba Loripes.orbiculatus
26.212753 25.931365 25.258931
Melinna.palmata Mytilaster.lineatus Oligochaeta
24.202974 21.932810 21.776133
Ampelisca.diadema Protodorvillea.kefersteini Salvatoria.clavata
20.058951 17.536359 17.003452
Polydora.ciliata Amphibalanus.improvisus Prionospio.cirrifera
15.689516 14.882244 13.192537
Sphaerosyllis.hystrix Rissoa.splendida Heteromastus.filiformis
12.713470 11.548802 10.221754
Lagis.koreni Microphthalmus.sp. Alitta.succinea
9.704792 9.395618 8.646146
Glycera.sp. Tricolia.pullus Platynereis.dumerilii
7.992056 7.955930 7.279560
Ampithoe.sp. Leiochone.leiopygos Syllis.gracilis
6.321173 6.218039 5.956444
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.1) <- names(top.sp.glms.lvm.zostera.1) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.1
Monocorophium acherusicum Cumella limicola Bittium reticulatum
38.328670 32.545296 31.182981
Apseudopsis ostroumovi Spio filicornis Chamelea gallina
30.854772 28.108956 26.230014
Capitella capitata Abra alba Loripes orbiculatus
26.212753 25.931365 25.258931
Melinna palmata Mytilaster lineatus Oligochaeta
24.202974 21.932810 21.776133
Ampelisca diadema Protodorvillea kefersteini Salvatoria clavata
20.058951 17.536359 17.003452
Polydora ciliata Amphibalanus improvisus Prionospio cirrifera
15.689516 14.882244 13.192537
Sphaerosyllis hystrix Rissoa splendida Heteromastus filiformis
12.713470 11.548802 10.221754
Lagis koreni Microphthalmus sp. Alitta succinea
9.704792 9.395618 8.646146
Glycera sp. Tricolia pullus Platynereis dumerilii
7.992056 7.955930 7.279560
Ampithoe sp. Leiochone leiopygos Syllis gracilis
6.321173 6.218039 5.956444
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.1 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.1)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.1))))
)
(plot.top.sp.glms.lvm.zostera.1 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.1,
mapping = aes(x = species, y = log_y_min(count), colour = station),
labs.legend = paste0("Z", as.numeric(unique(abnd.top.sp.glms.lvm.zostera.1$station))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Well this is a nightmarish plot, but more tolerable than the sand one - there are less species here, so at least it’s readable..
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.1 <- lapply(names(glms.lvm.zostera.1.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.1.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera2" etc.
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.1")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.1", "group_"))))
top.sp.abnd.glms.lvm.zostera.1 <- lapply(top.sp.abnd.glms.lvm.zostera.1, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.1 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station %in% c("Gradina", "Ropotamo") ~ 3)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.1.smry <- zoo.abnd.zostera.long.1 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.1 <- map2(top.sp.abnd.glms.lvm.zostera.1, zoo.abnd.zostera.long.1.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.1.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.1 <- map2(top.sp.abnd.glms.lvm.zostera.1, c(16, 4, 12), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
NA
In this case, the model shows which species exhibit a reaction based on the chosen groups - in other words, which species are more likely to be more/less abundant in each group.
I have to say, in the case of the seagrasses and case 1 clusters, there are much fewer species that exhibit a significant response - around 10 for each group.
The LRs are lower for groups 2 and 3 - not sure if this means anything, but for group 1 they are much much higher..
For group 1 (= Z1-Z2), the species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Oligochaeta, H. filiformis, Polydora ciliata, Prionospio cirrifera, R. membranacea, A. alba, A.diadema, M. acherusicum; and the only one with a significantly lower abundance - Chamelea gallina.
For group 2 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A.dadema. The species with lower abundance are: B. reticulatum, A. alba, Oligochaeta, S. clavata, P. ciliata, P. cirrifera, H. filiformis.
For group 3 (= Z4-Z5), the species with higher abundance are: Cumella limicola, Apseudopsis ostroumovi, Capitella capitata, Mytilaster lineatus, Loripes orbiculatus; less so, but still present - C. gallina, S. clavata. The species with lower abundance are: Abra alba, Melinna palmata (totally absent).
I’ll test each station as its own group, too (as I did before, with the classical multivariate methods) - I’m not sure how much I can trust this grouping (in particular group 3 is a bit far-fetched, if you ask me..).
LVM clusters - case 2 Check the model assumptions.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2, table = TRUE)
START SECTION 2
Plotting if overlay is TRUE
using grouping variable lvm.clusters.zostera.2 215 mean values were 0 and could
not be included in the log-plot
using grouping variable lvm.clusters.zostera.2 215 variance values were 0 and could not
be included in the log-plot
FINISHED SECTION 2
$mean
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 56.750 60.000 11.250 10.625 0.750
2 58.625 41.625 48.875 1.875 5.750
3 0.000 88.000 0.750 80.500 2.250
4 238.375 23.125 48.625 24.500 76.375
5 74.250 124.000 201.250 20.250 2.750
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata
1 8.625 12.50 46.750
2 27.000 38.25 6.375
3 2.250 0.00 0.000
4 23.000 0.00 0.250
5 5.250 20.00 11.250
Monocorophium.acherusicum Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi
1 3.500 5.375 0.125 0.000
2 4.375 9.500 0.750 0.125
3 78.250 9.500 0.000 0.000
4 1.000 10.125 9.875 3.875
5 2.250 4.000 39.750 48.250
Spio.filicornis Microdeutopus.gryllotalpa Abra.alba Cumella.limicola
1 1.250 3.625 6.125 0.00
2 1.625 1.625 6.500 0.75
3 38.250 4.250 0.000 0.00
4 1.250 1.875 0.750 10.25
5 0.000 14.000 1.750 6.50
Loripes.orbiculatus Parvicardium.exiguum Protodorvillea.kefersteini
1 0.625 2.625 0.250
2 1.625 3.000 1.250
3 0.500 5.500 0.000
4 8.750 2.875 9.875
5 4.500 1.500 1.250
Platynereis.dumerilii Nemertea Syllis.gracilis Alitta.succinea Amphibalanus.improvisus
1 1.875 1.000 0.000 7.125 5.500
2 1.250 3.750 3.250 0.000 1.125
3 9.250 0.500 0.000 0.000 0.250
4 1.625 3.125 0.875 0.000 0.250
5 3.750 1.500 7.250 1.000 0.500
Stenosoma.capito Lagis.koreni Schistomeringos.rudolphi Salvatoria.clavata
1 0.250 0.875 0.000 0.000
2 2.750 3.375 3.375 0.000
3 0.000 0.000 0.000 1.500
4 0.875 1.125 1.250 0.375
5 3.750 0.500 0.000 6.250
Leiochone.leiopygos Melinna.palmata Microphthalmus.sp. Kellia.suborbicularis
1 0.500 0.50 0.000 0.250
2 0.750 1.25 0.000 0.375
3 0.000 2.75 0.000 0.000
4 1.875 0.00 0.625 2.375
5 0.000 0.00 5.000 0.000
Nototropis.guttatus Chamelea.gallina Perinereis.cultrifera Sphaerosyllis.hystrix
1 0.625 0.000 0.750 0.000
2 0.875 0.125 0.125 0.375
3 0.000 0.000 0.000 0.000
4 0.125 2.000 0.625 1.375
5 2.500 1.250 2.000 1.250
Ampithoe.sp. Harmothoe.reticulata Phoronida Perioculodes.longimanus Rissoa.splendida
1 0.000 0.250 0.500 0.375 0.00
2 0.375 1.250 0.625 0.375 0.00
3 2.750 0.000 0.000 1.250 1.00
4 0.500 0.125 0.750 0.125 0.25
5 0.000 1.000 0.250 0.750 2.25
Diogenes.pugilator Iphinoe.tenella Platyhelminthes Exogone.naidina
1 0.25 0.125 1.125 0.000
2 0.50 1.250 0.250 0.125
3 0.75 0.000 0.000 2.250
4 0.25 0.250 0.250 0.250
5 0.75 0.000 0.000 0.000
Genetyllis.tuberculata Parthenina.interstincta Tritia.reticulata Cytharella.costulata
1 0.000 1.125 0.000 0.000
2 1.000 0.125 1.125 0.375
3 0.250 0.000 0.000 0.250
4 0.375 0.250 0.250 0.625
5 0.000 0.000 0.250 0.500
Syllis.hyalina Tricolia.pullus Turbellaria Harmothoe.imbricata Hydrobia.sp.
1 0.0 0.000 0.625 0.625 0.750
2 0.0 0.125 0.125 0.000 0.000
3 0.0 0.000 0.250 0.000 0.000
4 0.0 0.250 0.375 0.250 0.125
5 2.5 1.750 0.000 0.000 0.000
Lucinella.divaricata Nereis.pulsatoria Polititapes.aureus Upogebia.pusilla Glycera.sp.
1 0.000 0.875 0.000 0.00 0.375
2 0.000 0.000 0.750 0.00 0.375
3 0.000 0.000 0.000 0.00 0.000
4 0.875 0.000 0.125 0.00 0.000
5 0.000 0.000 0.000 1.75 0.000
Microphthalmus.fragilis Eunice.vittata Glycera.tridactyla Tritia.neritea
1 0.00 0.250 0.125 0.500
2 0.75 0.125 0.500 0.000
3 0.00 0.000 0.000 0.000
4 0.00 0.000 0.000 0.125
5 0.00 0.500 0.000 0.000
Chironomidae.larvae Hydrobia.acuta Micromaldane.ornithochaeta Cumella.pygmaea
1 0 0.5 0.000 0.000
2 0 0.0 0.000 0.000
3 0 0.0 0.250 0.000
4 0 0.0 0.375 0.375
5 1 0.0 0.000 0.000
Eurydice.dollfusi Hirudinea Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis
1 0.000 0.125 0.125 0.00 0.000
2 0.000 0.125 0.000 0.00 0.125
3 0.000 0.000 0.000 0.00 0.000
4 0.375 0.125 0.125 0.25 0.250
5 0.000 0.000 0.250 0.25 0.000
Brachynotus.sexdentatus Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa
1 0.125 0.00 0.125 0.000
2 0.000 0.00 0.000 0.125
3 0.000 0.00 0.000 0.000
4 0.125 0.25 0.000 0.125
5 0.000 0.00 0.250 0.000
Abra.sp. Actiniaria Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava
1 0.00 0.125 0.125 0.000 0.00
2 0.00 0.000 0.000 0.000 0.00
3 0.25 0.000 0.000 0.000 0.25
4 0.00 0.000 0.000 0.125 0.00
5 0.00 0.000 0.000 0.000 0.00
Gastrosaccus.sanctus Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis
1 0.000 0.125 0.000 0.000
2 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000
4 0.125 0.000 0.125 0.125
5 0.000 0.000 0.000 0.000
Maldane.glebifex Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.125 0.125 0.125 0.000
2 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000
4 0.000 0.000 0.000 0.125
5 0.000 0.000 0.000 0.000
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.125 0.00 0.00 0.00
2 0.000 0.00 0.00 0.00
3 0.000 0.25 0.25 0.25
4 0.000 0.00 0.00 0.00
5 0.000 0.00 0.00 0.00
Thracia.phaseolina
1 0.000
2 0.000
3 0.000
4 0.125
5 0.000
$var
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 4684.5000 4366.8571 49.07143 49.410714 1.928571
2 975.4107 1652.8393 643.83929 4.410714 58.785714
3 0.0000 3094.0000 2.25000 395.000000 1.583333
4 25180.5536 452.6964 3778.26786 735.142857 877.125000
5 3278.2500 3103.3333 18570.91667 92.916667 2.916667
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata
1 77.125000 316.5714 1167.6428571
2 110.285714 1275.0714 44.2678571
3 1.583333 0.0000 0.0000000
4 136.000000 0.0000 0.2142857
5 32.916667 183.3333 34.2500000
Monocorophium.acherusicum Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi
1 22.000000 12.839286 0.125000 0.000000
2 15.982143 559.714286 3.071429 0.125000
3 3980.916667 11.666667 0.000000 0.000000
4 1.428571 95.553571 82.982143 9.839286
5 4.250000 6.666667 566.250000 196.250000
Spio.filicornis Microdeutopus.gryllotalpa Abra.alba Cumella.limicola
1 3.357143 19.410714 13.267857 0.000000
2 1.982143 1.410714 38.857143 2.214286
3 306.250000 16.250000 0.000000 0.000000
4 3.357143 2.696429 1.071429 30.500000
5 0.000000 40.666667 1.583333 23.000000
Loripes.orbiculatus Parvicardium.exiguum Protodorvillea.kefersteini
1 0.8392857 7.125000 0.50000
2 3.6964286 33.142857 4.50000
3 1.0000000 3.666667 0.00000
4 16.5000000 8.410714 16.69643
5 7.0000000 1.666667 2.25000
Platynereis.dumerilii Nemertea Syllis.gracilis Alitta.succinea
1 12.125000 1.7142857 0.000000 35.5535714
2 2.500000 20.2142857 9.928571 0.0000000
3 32.916667 0.3333333 0.000000 0.0000000
4 2.553571 3.8392857 2.982143 0.0000000
5 10.916667 3.6666667 22.916667 0.6666667
Amphibalanus.improvisus Stenosoma.capito Lagis.koreni Schistomeringos.rudolphi
1 38.0000000 0.2142857 1.5535714 0.000000
2 1.2678571 38.5000000 8.8392857 17.696429
3 0.2500000 0.0000000 0.0000000 0.000000
4 0.2142857 0.9821429 0.9821429 3.642857
5 1.0000000 12.9166667 0.3333333 0.000000
Salvatoria.clavata Leiochone.leiopygos Melinna.palmata Microphthalmus.sp.
1 0.0000000 0.5714286 0.5714286 0.000
2 0.0000000 0.7857143 1.6428571 0.000
3 0.3333333 0.0000000 0.2500000 0.000
4 1.1250000 3.8392857 0.0000000 3.125
5 14.9166667 0.0000000 0.0000000 22.000
Kellia.suborbicularis Nototropis.guttatus Chamelea.gallina Perinereis.cultrifera
1 0.2142857 1.982143 0.000000 4.5000000
2 1.1250000 1.553571 0.125000 0.1250000
3 0.0000000 0.000000 0.000000 0.0000000
4 13.6964286 0.125000 1.714286 0.8392857
5 0.0000000 14.333333 1.583333 3.3333333
Sphaerosyllis.hystrix Ampithoe.sp. Harmothoe.reticulata Phoronida
1 0.0000000 0.0000000 0.2142857 0.5714286
2 0.5535714 0.5535714 0.7857143 0.5535714
3 0.0000000 11.5833333 0.0000000 0.0000000
4 2.2678571 0.5714286 0.1250000 1.0714286
5 0.9166667 0.0000000 2.0000000 0.2500000
Perioculodes.longimanus Rissoa.splendida Diogenes.pugilator Iphinoe.tenella
1 0.5535714 0.0000000 0.2142857 0.125000
2 0.5535714 0.0000000 0.5714286 3.928571
3 3.5833333 2.0000000 0.2500000 0.000000
4 0.1250000 0.2142857 0.2142857 0.500000
5 0.9166667 6.9166667 0.9166667 0.000000
Platyhelminthes Exogone.naidina Genetyllis.tuberculata Parthenina.interstincta
1 2.1250000 0.000 0.0000000 4.982143
2 0.5000000 0.125 4.2857143 0.125000
3 0.0000000 8.250 0.2500000 0.000000
4 0.2142857 0.500 0.5535714 0.500000
5 0.0000000 0.000 0.0000000 0.000000
Tritia.reticulata Cytharella.costulata Syllis.hyalina Tricolia.pullus Turbellaria
1 0.0000000 0.0000000 0 0.0000000 1.982143
2 2.4107143 0.2678571 0 0.1250000 0.125000
3 0.0000000 0.2500000 0 0.0000000 0.250000
4 0.2142857 0.8392857 0 0.2142857 1.125000
5 0.2500000 0.3333333 25 2.9166667 0.000000
Harmothoe.imbricata Hydrobia.sp. Lucinella.divaricata Nereis.pulsatoria
1 1.1250000 2.214286 0.000000 2.696429
2 0.0000000 0.000000 0.000000 0.000000
3 0.0000000 0.000000 0.000000 0.000000
4 0.2142857 0.125000 3.267857 0.000000
5 0.0000000 0.000000 0.000000 0.000000
Polititapes.aureus Upogebia.pusilla Glycera.sp. Microphthalmus.fragilis Eunice.vittata
1 0.000000 0.000000 0.5535714 0.0 0.500
2 1.071429 0.000000 0.2678571 4.5 0.125
3 0.000000 0.000000 0.0000000 0.0 0.000
4 0.125000 0.000000 0.0000000 0.0 0.000
5 0.000000 5.583333 0.0000000 0.0 1.000
Glycera.tridactyla Tritia.neritea Chironomidae.larvae Hydrobia.acuta
1 0.125000 1.142857 0 2
2 1.142857 0.000000 0 0
3 0.000000 0.000000 0 0
4 0.000000 0.125000 0 0
5 0.000000 0.000000 2 0
Micromaldane.ornithochaeta Cumella.pygmaea Eurydice.dollfusi Hirudinea
1 0.0000000 0.0000000 0.0000000 0.125
2 0.0000000 0.0000000 0.0000000 0.125
3 0.2500000 0.0000000 0.0000000 0.000
4 0.5535714 0.5535714 0.5535714 0.125
5 0.0000000 0.0000000 0.0000000 0.000
Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis Brachynotus.sexdentatus
1 0.125 0.0000000 0.0000000 0.125
2 0.000 0.0000000 0.1250000 0.000
3 0.000 0.0000000 0.0000000 0.000
4 0.125 0.2142857 0.2142857 0.125
5 0.250 0.2500000 0.0000000 0.000
Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa Abra.sp. Actiniaria
1 0.0000000 0.125 0.000 0.00 0.125
2 0.0000000 0.000 0.125 0.00 0.000
3 0.0000000 0.000 0.000 0.25 0.000
4 0.2142857 0.000 0.125 0.00 0.000
5 0.0000000 0.250 0.000 0.00 0.000
Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava Gastrosaccus.sanctus
1 0.125 0.000 0.00 0.000
2 0.000 0.000 0.00 0.000
3 0.000 0.000 0.25 0.000
4 0.000 0.125 0.00 0.125
5 0.000 0.000 0.00 0.000
Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis Maldane.glebifex
1 0.125 0.000 0.000 0.125
2 0.000 0.000 0.000 0.000
3 0.000 0.000 0.000 0.000
4 0.000 0.125 0.125 0.000
5 0.000 0.000 0.000 0.000
Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.125 0.125 0.000
2 0.000 0.000 0.000
3 0.000 0.000 0.000
4 0.000 0.000 0.125
5 0.000 0.000 0.000
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.125 0.00 0.00 0.00
2 0.000 0.00 0.00 0.00
3 0.000 0.25 0.25 0.25
4 0.000 0.00 0.00 0.00
5 0.000 0.00 0.00 0.00
Thracia.phaseolina
1 0.000
2 0.000
3 0.000
4 0.125
5 0.000
It’s not perfect, but it’s not too terrible either. I think it’s a little worse than the case 1 fit.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.2 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.2,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.2)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.2, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(2,2))
# lapply(2:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(2, 2))
Save the model!
write_rds(glms.lvm.zostera.2,
here(save.dir, "glms_lvm_zostera_2.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.2.summary <- summary(glms.lvm.zostera.2,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
Resampling run 0 finished. Time elapsed: 0.00 min ...
Resampling run 100 finished. Time elapsed: 0.42 min ...
Resampling run 200 finished. Time elapsed: 0.84 min ...
Resampling run 300 finished. Time elapsed: 1.26 min ...
Resampling run 400 finished. Time elapsed: 1.67 min ...
Resampling run 500 finished. Time elapsed: 2.09 min ...
Resampling run 600 finished. Time elapsed: 2.51 min ...
Resampling run 700 finished. Time elapsed: 2.93 min ...
Resampling run 800 finished. Time elapsed: 3.35 min ...
Resampling run 900 finished. Time elapsed: 3.77 min ...
Time elapsed: 0 hr 4 min 11 sec
Test statistics:
LR value Pr(>LR)
(Intercept) 1080.0 0.001 ***
lvm.clusters.zostera.22 256.8 0.001 ***
lvm.clusters.zostera.23 310.9 0.001 ***
lvm.clusters.zostera.24 471.9 0.001 ***
lvm.clusters.zostera.25 346.6 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate test statistic:
(Intercept) lvm.clusters.zostera.22
LR value Pr(>LR) LR value Pr(>LR)
Abra.alba 31.096 0.001 0.023 1.000
Abra.sp. 4.677 0.714 0.000 1.000
Actiniaria 2.341 0.978 1.386 0.990
Alitta.succinea 25.658 0.001 23.057 0.002
Ampelisca.diadema 54.649 0.001 10.416 0.072
Amphibalanus.improvisus 23.546 0.001 8.177 0.181
Ampithoe.sp. 8.217 0.202 3.482 0.823
Anadara.kagoshimensis 2.341 0.978 1.386 0.990
Apherusa.bispinosa 3.770 0.854 0.000 1.000
Apseudopsis.ostroumovi 14.963 0.013 1.380 0.992
Bittium.reticulatum 126.101 0.001 0.004 1.000
Brachynotus.sexdentatus 3.387 0.910 1.386 0.990
Capitella.capitata 5.997 0.502 2.862 0.912
Capitella.minima 132.041 0.001 0.433 1.000
Chamelea.gallina 15.617 0.009 1.386 0.990
Chironomidae.larvae 6.560 0.425 0.000 1.000
Cumella.limicola 13.674 0.027 7.694 0.208
Cumella.pygmaea 5.165 0.620 0.000 1.000
Cytharella.costulata 12.759 0.036 4.159 0.675
Diogenes.pugilator 6.380 0.439 0.680 0.999
Eteone.flava 4.677 0.714 0.000 1.000
Eunice.vittata 1.521 0.992 0.208 1.000
Eurydice.dollfusi 5.165 0.620 0.000 1.000
Exogone.naidina 5.874 0.506 1.250 0.995
Gastrosaccus.sanctus 3.770 0.854 0.000 1.000
Genetyllis.tuberculata 6.742 0.398 5.746 0.444
Glycera.sp. 2.965 0.944 0.000 1.000
Glycera.tridactyla 2.776 0.955 1.033 0.998
Glycera.unicornis 2.341 0.978 1.386 0.990
Harmothoe.imbricata 0.759 1.000 5.561 0.476
Harmothoe.reticulata 6.385 0.439 5.794 0.444
Heteromastus.filiformis 62.249 0.001 10.679 0.060
Hirudinea 4.137 0.812 0.000 1.000
Hydrobia.acuta 0.171 1.000 1.494 0.983
Hydrobia.sp. 0.116 1.000 3.709 0.788
Iphinoe.tenella 4.775 0.705 4.354 0.646
Kellia.suborbicularis 1.985 0.978 0.104 1.000
Lagis.koreni 0.095 1.000 5.999 0.421
Leiochone.leiopygos 1.899 0.978 0.336 1.000
Lentidium.mediterraneum 6.813 0.388 0.000 1.000
Lepidochitona.cinerea 3.770 0.854 0.000 1.000
Loripes.orbiculatus 1.160 0.995 3.172 0.869
Lucinella.divaricata 3.682 0.861 0.000 1.000
Magelona.papillicornis 3.770 0.854 0.000 1.000
Maldane.glebifex 2.341 0.978 1.386 0.990
Melinna.palmata 2.455 0.978 2.657 0.931
Microdeutopus.gryllotalpa 14.797 0.013 2.603 0.933
Micromaldane.ornithochaeta 6.542 0.426 0.000 1.000
Micronephthys.stammeri 3.925 0.837 1.386 0.990
Microphthalmus.fragilis 1.716 0.983 1.494 0.983
Microphthalmus.sp. 5.680 0.536 0.000 1.000
Monocorophium.acherusicum 12.799 0.036 0.213 1.000
Mytilaster.lineatus 0.385 1.000 12.953 0.032
Mytilus.galloprovincialis 2.341 0.978 1.386 0.990
Nemertea 0.000 1.000 5.479 0.494
Nephtys.cirrosa 5.423 0.582 1.386 0.990
Nephtys.kersivalensis 2.341 0.978 1.386 0.990
Nereis.perivisceralis 3.770 0.854 0.000 1.000
Nereis.pulsatoria 0.020 1.000 3.251 0.853
Nototropis.guttatus 0.402 1.000 0.116 1.000
Oligochaeta 58.093 0.001 9.386 0.110
Paradoneis.harpagonea 2.341 0.978 1.386 0.990
Parthenina.interstincta 0.017 1.000 1.662 0.983
Parvicardium.exiguum 6.468 0.429 0.063 1.000
Perinereis.cultrifera 0.177 1.000 2.035 0.967
Perioculodes.longimanus 2.180 0.978 0.000 1.000
Phoronida 2.237 0.978 0.110 1.000
Phyllodoce.sp. 4.677 0.714 0.000 1.000
Platyhelminthes 0.058 1.000 2.602 0.933
Platynereis.dumerilii 1.999 0.978 0.374 1.000
Polititapes.aureus 9.599 0.112 6.933 0.288
Polychaeta.larvae 4.677 0.714 0.000 1.000
Polydora.ciliata 91.477 0.001 12.416 0.033
Polygordius.neapolitanus 4.677 0.714 0.000 1.000
Prionospio.cirrifera 49.773 0.001 4.922 0.568
Protodorvillea.kefersteini 5.371 0.585 4.879 0.568
Pseudocuma.longicorne 4.585 0.729 1.386 0.990
Rissoa.membranacea 20.908 0.001 1.019 0.998
Rissoa.splendida 10.372 0.071 0.000 1.000
Salvatoria.clavata 14.727 0.013 0.000 1.000
Schistomeringos.rudolphi 6.673 0.412 8.799 0.141
Sphaerosyllis.hystrix 12.303 0.039 4.024 0.720
Spio.filicornis 0.320 1.000 0.246 1.000
Spisula.subtruncata 8.280 0.194 0.000 1.000
Stenosoma.capito 3.207 0.929 6.661 0.330
Syllis.gracilis 11.730 0.042 17.831 0.005
Syllis.hyalina 2.396 0.978 0.000 1.000
Tellina.tenuis 7.750 0.294 1.386 0.990
Thracia.phaseolina 3.770 0.854 0.000 1.000
Tricolia.pullus 11.452 0.047 1.380 0.992
Tritia.neritea 0.707 1.000 3.479 0.823
Tritia.reticulata 9.399 0.132 8.301 0.174
Turbellaria 0.265 1.000 1.125 0.998
Upogebia.pusilla 5.886 0.505 0.000 1.000
lvm.clusters.zostera.23 lvm.clusters.zostera.24
LR value Pr(>LR) LR value
Abra.alba 17.297 0.003 13.450
Abra.sp. 2.197 0.783 0.000
Actiniaria 0.811 0.947 1.386
Alitta.succinea 15.046 0.008 23.057
Ampelisca.diadema 12.759 0.014 3.426
Amphibalanus.improvisus 10.584 0.036 16.614
Ampithoe.sp. 9.419 0.062 4.372
Anadara.kagoshimensis 0.811 0.947 1.386
Apherusa.bispinosa 0.000 1.000 1.386
Apseudopsis.ostroumovi 0.000 0.969 26.746
Bittium.reticulatum 24.917 0.001 6.987
Brachynotus.sexdentatus 0.811 0.947 0.000
Capitella.capitata 0.773 0.950 17.684
Capitella.minima 0.333 0.969 2.753
Chamelea.gallina 0.000 1.000 19.316
Chironomidae.larvae 0.000 0.969 0.000
Cumella.limicola 0.000 1.000 31.096
Cumella.pygmaea 0.000 1.000 3.167
Cytharella.costulata 2.197 0.783 6.931
Diogenes.pugilator 1.483 0.898 0.000
Eteone.flava 2.197 0.783 0.000
Eunice.vittata 1.199 0.898 2.081
Eurydice.dollfusi 0.000 1.000 3.167
Exogone.naidina 6.160 0.228 2.262
Gastrosaccus.sanctus 0.000 1.000 1.386
Genetyllis.tuberculata 1.876 0.857 3.160
Glycera.sp. 2.428 0.738 4.124
Glycera.tridactyla 0.718 0.958 1.242
Glycera.unicornis 0.811 0.947 1.386
Harmothoe.imbricata 3.293 0.613 0.994
Harmothoe.reticulata 1.622 0.893 0.340
Heteromastus.filiformis 7.032 0.170 8.276
Hirudinea 0.811 0.947 0.000
Hydrobia.acuta 0.884 0.941 1.494
Hydrobia.sp. 2.194 0.783 1.407
Iphinoe.tenella 0.762 0.952 0.276
Kellia.suborbicularis 1.257 0.898 3.343
Lagis.koreni 4.950 0.375 0.180
Leiochone.leiopygos 3.040 0.670 4.486
Lentidium.mediterraneum 0.000 1.000 2.773
Lepidochitona.cinerea 0.000 1.000 1.386
Loripes.orbiculatus 0.068 0.969 23.773
Lucinella.divaricata 0.000 1.000 3.204
Magelona.papillicornis 0.000 1.000 1.386
Maldane.glebifex 0.811 0.947 1.386
Melinna.palmata 9.495 0.062 5.545
Microdeutopus.gryllotalpa 0.089 0.969 1.848
Micromaldane.ornithochaeta 2.067 0.809 3.568
Micronephthys.stammeri 0.811 0.947 1.386
Microphthalmus.fragilis 0.000 0.969 0.000
Microphthalmus.sp. 0.000 0.969 3.972
Monocorophium.acherusicum 20.414 0.002 4.667
Mytilaster.lineatus 2.713 0.724 35.697
Mytilus.galloprovincialis 0.811 0.947 1.386
Nemertea 0.608 0.958 4.090
Nephtys.cirrosa 0.000 1.000 1.386
Nephtys.kersivalensis 0.811 0.947 1.386
Nereis.perivisceralis 0.000 1.000 1.386
Nereis.pulsatoria 1.949 0.850 3.251
Nototropis.guttatus 2.528 0.725 1.537
Oligochaeta 11.422 0.027 9.333
Paradoneis.harpagonea 0.811 0.947 1.386
Parthenina.interstincta 2.228 0.775 0.968
Parvicardium.exiguum 1.435 0.898 0.029
Perinereis.cultrifera 2.924 0.685 0.035
Perioculodes.longimanus 1.591 0.898 0.850
Phoronida 3.210 0.631 0.396
Phyllodoce.sp. 2.197 0.783 0.000
Platyhelminthes 4.577 0.403 2.602
Platynereis.dumerilii 4.688 0.403 0.049
Polititapes.aureus 0.000 1.000 1.367
Polychaeta.larvae 2.197 0.783 0.000
Polydora.ciliata 24.008 0.001 30.614
Polygordius.neapolitanus 2.197 0.783 0.000
Prionospio.cirrifera 18.028 0.002 29.724
Protodorvillea.kefersteini 1.584 0.898 20.929
Pseudocuma.longicorne 0.811 0.947 0.000
Rissoa.membranacea 0.725 0.958 1.256
Rissoa.splendida 7.113 0.167 2.626
Salvatoria.clavata 12.384 0.017 4.104
Schistomeringos.rudolphi 0.000 1.000 6.233
Sphaerosyllis.hystrix 0.000 0.969 11.891
Spio.filicornis 19.150 0.002 0.000
Spisula.subtruncata 0.000 0.975 2.773
Stenosoma.capito 1.424 0.898 1.680
Syllis.gracilis 0.000 1.000 8.039
Syllis.hyalina 0.000 0.969 0.000
Tellina.tenuis 0.000 1.000 2.773
Thracia.phaseolina 0.000 1.000 1.386
Tricolia.pullus 0.000 0.969 2.745
Tritia.neritea 2.051 0.812 1.033
Tritia.reticulata 0.000 0.969 2.586
Turbellaria 0.285 0.969 0.158
Upogebia.pusilla 0.000 1.000 0.000
lvm.clusters.zostera.25
Pr(>LR) LR value Pr(>LR)
Abra.alba 0.018 4.346 0.669
Abra.sp. 1.000 0.000 1.000
Actiniaria 0.996 0.811 0.999
Alitta.succinea 0.002 6.818 0.295
Ampelisca.diadema 0.806 1.516 0.991
Amphibalanus.improvisus 0.005 8.450 0.159
Ampithoe.sp. 0.642 0.000 1.000
Anadara.kagoshimensis 0.996 0.811 0.999
Apherusa.bispinosa 0.995 0.000 0.999
Apseudopsis.ostroumovi 0.001 39.178 0.001
Bittium.reticulatum 0.249 0.201 0.999
Brachynotus.sexdentatus 1.000 0.811 0.999
Capitella.capitata 0.005 22.609 0.001
Capitella.minima 0.934 1.211 0.994
Chamelea.gallina 0.005 10.986 0.050
Chironomidae.larvae 1.000 5.456 0.476
Cumella.limicola 0.001 26.312 0.001
Cumella.pygmaea 0.849 0.000 1.000
Cytharella.costulata 0.256 4.394 0.664
Diogenes.pugilator 1.000 1.483 0.991
Eteone.flava 1.000 0.000 1.000
Eunice.vittata 0.983 0.210 0.999
Eurydice.dollfusi 0.849 0.000 0.999
Exogone.naidina 0.972 0.000 0.999
Gastrosaccus.sanctus 0.996 0.000 1.000
Genetyllis.tuberculata 0.850 0.000 1.000
Glycera.sp. 0.721 2.428 0.938
Glycera.tridactyla 0.997 0.718 0.999
Glycera.unicornis 0.996 0.811 0.999
Harmothoe.imbricata 0.998 3.293 0.805
Harmothoe.reticulata 1.000 2.773 0.880
Heteromastus.filiformis 0.172 1.330 0.994
Hirudinea 1.000 0.811 0.999
Hydrobia.acuta 0.993 0.884 0.999
Hydrobia.sp. 0.995 2.194 0.953
Iphinoe.tenella 1.000 0.762 0.999
Kellia.suborbicularis 0.818 1.257 0.994
Lagis.koreni 1.000 0.420 0.999
Leiochone.leiopygos 0.624 3.040 0.837
Lentidium.mediterraneum 0.921 0.000 1.000
Lepidochitona.cinerea 0.995 0.000 1.000
Loripes.orbiculatus 0.002 12.893 0.023
Lucinella.divaricata 0.844 0.000 1.000
Magelona.papillicornis 0.996 0.000 1.000
Maldane.glebifex 0.996 0.811 0.999
Melinna.palmata 0.470 3.244 0.807
Microdeutopus.gryllotalpa 0.989 6.470 0.331
Micromaldane.ornithochaeta 0.786 0.000 1.000
Micronephthys.stammeri 0.995 0.236 0.999
Microphthalmus.fragilis 1.000 0.000 1.000
Microphthalmus.sp. 0.731 7.235 0.253
Monocorophium.acherusicum 0.605 0.479 0.999
Mytilaster.lineatus 0.001 3.994 0.738
Mytilus.galloprovincialis 0.996 0.811 0.999
Nemertea 0.722 0.324 0.999
Nephtys.cirrosa 0.995 0.000 0.999
Nephtys.kersivalensis 0.996 0.811 0.999
Nereis.perivisceralis 0.996 0.000 1.000
Nereis.pulsatoria 0.831 1.949 0.972
Nototropis.guttatus 0.992 1.502 0.991
Oligochaeta 0.100 21.743 0.001
Paradoneis.harpagonea 0.996 0.811 0.999
Parthenina.interstincta 0.998 2.228 0.951
Parvicardium.exiguum 1.000 0.610 0.999
Perinereis.cultrifera 1.000 0.836 0.999
Perioculodes.longimanus 0.998 0.467 0.999
Phoronida 1.000 0.433 0.999
Phyllodoce.sp. 1.000 0.000 1.000
Platyhelminthes 0.954 4.577 0.632
Platynereis.dumerilii 1.000 0.899 0.999
Polititapes.aureus 0.996 0.000 1.000
Polychaeta.larvae 1.000 0.000 1.000
Polydora.ciliata 0.001 5.015 0.546
Polygordius.neapolitanus 1.000 0.000 1.000
Prionospio.cirrifera 0.001 0.688 0.999
Protodorvillea.kefersteini 0.004 3.630 0.780
Pseudocuma.longicorne 1.000 0.236 0.999
Rissoa.membranacea 0.997 0.171 0.999
Rissoa.splendida 0.950 11.194 0.047
Salvatoria.clavata 0.722 23.844 0.001
Schistomeringos.rudolphi 0.366 0.000 0.999
Sphaerosyllis.hystrix 0.026 9.470 0.103
Spio.filicornis 1.000 6.560 0.319
Spisula.subtruncata 0.921 2.197 0.953
Stenosoma.capito 0.989 7.067 0.269
Syllis.gracilis 0.175 21.341 0.001
Syllis.hyalina 1.000 2.507 0.929
Tellina.tenuis 0.931 0.000 1.000
Thracia.phaseolina 0.995 0.000 0.999
Tricolia.pullus 0.934 11.106 0.048
Tritia.neritea 0.998 2.051 0.963
Tritia.reticulata 0.954 2.073 0.961
Turbellaria 1.000 2.003 0.967
Upogebia.pusilla 1.000 5.769 0.409
Arguments: with 999 resampling iterations using pit.trap resampling and response assumed to be uncorrelated
Likelihood Ratio statistic: 1223, p-value: 0.001
Univariate test statistic:
Abra.alba Abra.sp. Actiniaria Alitta.succinea Ampelisca.diadema
LR value 27.461 4.159 2.773 36.585 30.570
Pr(>LR) 0.002 0.961 0.981 0.001 0.001
Amphibalanus.improvisus Ampithoe.sp. Anadara.kagoshimensis Apherusa.bispinosa
LR value 23.408 12.209 2.773 2.773
Pr(>LR) 0.008 0.282 0.981 0.981
Apseudopsis.ostroumovi Bittium.reticulatum Brachynotus.sexdentatus
LR value 52.794 34.193 2.773
Pr(>LR) 0.001 0.001 0.981
Capitella.capitata Capitella.minima Chamelea.gallina Chironomidae.larvae
LR value 33.199 7.227 28.525 9.651
Pr(>LR) 0.001 0.781 0.001 0.513
Cumella.limicola Cumella.pygmaea Cytharella.costulata Diogenes.pugilator
LR value 41.328 6.118 7.354 3.059
Pr(>LR) 0.001 0.840 0.775 0.981
Eteone.flava Eunice.vittata Eurydice.dollfusi Exogone.naidina
LR value 4.159 3.779 6.118 8.477
Pr(>LR) 0.961 0.967 0.840 0.642
Gastrosaccus.sanctus Genetyllis.tuberculata Glycera.sp. Glycera.tridactyla
LR value 2.773 7.712 7.992 5.516
Pr(>LR) 0.981 0.739 0.704 0.903
Glycera.unicornis Harmothoe.imbricata Harmothoe.reticulata
LR value 2.773 8.753 15.210
Pr(>LR) 0.981 0.597 0.093
Heteromastus.filiformis Hirudinea Hydrobia.acuta Hydrobia.sp. Iphinoe.tenella
LR value 28.142 1.726 2.932 5.959 9.824
Pr(>LR) 0.002 0.981 0.981 0.860 0.484
Kellia.suborbicularis Lagis.koreni Leiochone.leiopygos Lentidium.mediterraneum
LR value 9.077 16.471 15.331 5.545
Pr(>LR) 0.579 0.076 0.091 0.897
Lepidochitona.cinerea Loripes.orbiculatus Lucinella.divaricata
LR value 2.773 31.348 6.172
Pr(>LR) 0.981 0.001 0.830
Magelona.papillicornis Maldane.glebifex Melinna.palmata
LR value 2.773 2.773 26.860
Pr(>LR) 0.981 0.981 0.002
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta Micronephthys.stammeri
LR value 16.350 6.655 4.159
Pr(>LR) 0.076 0.812 0.961
Microphthalmus.fragilis Microphthalmus.sp. Monocorophium.acherusicum
LR value 2.931 11.855 39.891
Pr(>LR) 0.981 0.306 0.001
Mytilaster.lineatus Mytilus.galloprovincialis Nemertea Nephtys.cirrosa
LR value 47.863 2.773 10.046 2.773
Pr(>LR) 0.001 0.981 0.484 0.981
Nephtys.kersivalensis Nereis.perivisceralis Nereis.pulsatoria
LR value 2.773 2.773 6.238
Pr(>LR) 0.981 0.981 0.823
Nototropis.guttatus Oligochaeta Paradoneis.harpagonea Parthenina.interstincta
LR value 7.210 36.765 2.773 4.546
Pr(>LR) 0.781 0.001 0.981 0.951
Parvicardium.exiguum Perinereis.cultrifera Perioculodes.longimanus Phoronida
LR value 2.973 6.689 4.365 5.483
Pr(>LR) 0.981 0.812 0.961 0.906
Phyllodoce.sp. Platyhelminthes Platynereis.dumerilii Polititapes.aureus
LR value 4.159 8.508 8.910 10.998
Pr(>LR) 0.961 0.642 0.597 0.372
Polychaeta.larvae Polydora.ciliata Polygordius.neapolitanus
LR value 4.159 41.905 4.159
Pr(>LR) 0.961 0.001 0.961
Prionospio.cirrifera Protodorvillea.kefersteini Pseudocuma.longicorne
LR value 46.273 30.140 3.112
Pr(>LR) 0.001 0.001 0.981
Rissoa.membranacea Rissoa.splendida Salvatoria.clavata Schistomeringos.rudolphi
LR value 2.720 17.118 32.959 13.822
Pr(>LR) 0.981 0.053 0.001 0.172
Sphaerosyllis.hystrix Spio.filicornis Spisula.subtruncata Stenosoma.capito
LR value 16.760 34.879 5.885 13.850
Pr(>LR) 0.068 0.001 0.863 0.172
Syllis.gracilis Syllis.hyalina Tellina.tenuis Thracia.phaseolina
LR value 28.543 4.555 4.499 2.773
Pr(>LR) 0.001 0.951 0.951 0.981
Tricolia.pullus Tritia.neritea Tritia.reticulata Turbellaria Upogebia.pusilla
LR value 14.668 5.516 11.233 2.652 10.026
Pr(>LR) 0.118 0.903 0.353 0.981 0.484
Arguments:
Test statistics calculated assuming response assumed to be uncorrelated
P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).
The factor is highly significant according to the models.
Again, save the summary for safekeeping, but also run an anova.
write_rds(glms.lvm.zostera.2.summary,
here(save.dir, "glms_lvm_zostera_2_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.2.aov <- anova.manyglm(glms.lvm.zostera.2,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.24 minutes...
Resampling run 300 finished. Time elapsed: 0.36 minutes...
Resampling run 400 finished. Time elapsed: 0.47 minutes...
Resampling run 500 finished. Time elapsed: 0.59 minutes...
Resampling run 600 finished. Time elapsed: 0.71 minutes...
Resampling run 700 finished. Time elapsed: 0.83 minutes...
Resampling run 800 finished. Time elapsed: 0.96 minutes...
Resampling run 900 finished. Time elapsed: 1.09 minutes...
Time elapsed: 0 hr 1 min 12 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.2,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.2 27 4 1223 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 27.461 0.003 4.159 1.000 2.773 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 36.585 0.001 30.57 0.001
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 23.408 0.009 12.209 0.430
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 52.794 0.001 34.193 0.001
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 33.199 0.001
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.227 0.944 28.525 0.002
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 9.651 0.743 41.328 0.001
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 6.118 0.976 7.354 0.936
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 3.059 1.000 4.159 1.000 3.779
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.000 6.118 0.976 8.477 0.887
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 7.712 0.927
Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.992 0.921 5.516 0.996
Glycera.unicornis Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 8.753 0.855
Harmothoe.reticulata Heteromastus.filiformis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 15.21 0.185 28.142 0.002
Hirudinea Hydrobia.acuta Hydrobia.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.726 1.000 2.932 1.000 5.959 0.981
Iphinoe.tenella Kellia.suborbicularis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 9.824 0.727 9.077 0.837
Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 16.471 0.135 15.331 0.183
Lentidium.mediterraneum Lepidochitona.cinerea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 5.545 0.993 2.773 1.000
Loripes.orbiculatus Lucinella.divaricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 31.348 0.001 6.172 0.976
Magelona.papillicornis Maldane.glebifex
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Melinna.palmata Microdeutopus.gryllotalpa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 26.86 0.003 16.35 0.135
Micromaldane.ornithochaeta Micronephthys.stammeri
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 6.655 0.964 4.159
Microphthalmus.fragilis Microphthalmus.sp.
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 1.000 2.931 1.000 11.855
Monocorophium.acherusicum Mytilaster.lineatus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 0.463 39.891 0.001 47.863
Mytilus.galloprovincialis Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 0.001 2.773 1.000 10.046 0.715
Nephtys.cirrosa Nephtys.kersivalensis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 2.773 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 6.238 0.973
Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 7.21 0.944 36.765 0.001
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.773 1.000 4.546 0.998
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 2.973 1.000 6.689 0.964
Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.365 0.999 5.483 0.996
Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 1.000 8.508 0.887
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 8.91 0.850 10.998 0.573
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 1.000 41.905 0.001
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 4.159 1.000 46.273 0.001
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 30.14 0.002 3.112
Rissoa.membranacea Rissoa.splendida
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.000 2.72 1.000 17.118 0.115
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 32.959 0.001 13.822 0.294
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 16.76 0.128 34.879 0.001
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 5.885 0.981 13.85 0.294
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 28.543 0.002 4.555 0.998 4.499
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 0.998 2.773 1.000 14.668 0.221
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.2 5.516 0.996 11.233 0.554 2.652
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.2 1.000 10.026 0.715
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Save the ANOVA, too.
write_rds(glms.lvm.zostera.2.aov,
here(save.dir, "glms_lvm_zostera_2_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.2 <- top_n_sp_glm(glms.lvm.zostera.2.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.799"
Apseudopsis.ostroumovi Mytilaster.lineatus Prionospio.cirrifera
52.793744 47.863335 46.272692
Polydora.ciliata Cumella.limicola Monocorophium.acherusicum
41.905219 41.327723 39.891356
Oligochaeta Alitta.succinea Spio.filicornis
36.764590 36.584687 34.878811
Bittium.reticulatum Capitella.capitata Salvatoria.clavata
34.192799 33.198554 32.958829
Loripes.orbiculatus Ampelisca.diadema Protodorvillea.kefersteini
31.347733 30.569510 30.139960
Syllis.gracilis Chamelea.gallina Heteromastus.filiformis
28.542915 28.524575 28.142266
Abra.alba Melinna.palmata Amphibalanus.improvisus
27.461185 26.859546 23.408099
Rissoa.splendida Sphaerosyllis.hystrix Lagis.koreni
17.117581 16.760313 16.471224
Microdeutopus.gryllotalpa Leiochone.leiopygos Harmothoe.reticulata
16.349670 15.331303 15.209809
Tricolia.pullus Stenosoma.capito Schistomeringos.rudolphi
14.667741 13.849706 13.822279
Ampithoe.sp. Microphthalmus.sp. Tritia.reticulata
12.208710 11.854688 11.232600
Polititapes.aureus Nemertea Upogebia.pusilla
10.998455 10.045834 10.025894
Iphinoe.tenella Chironomidae.larvae Kellia.suborbicularis
9.824385 9.650905 9.077281
Platynereis.dumerilii
8.910425
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.2) <- names(top.sp.glms.lvm.zostera.2) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.2
Apseudopsis ostroumovi Mytilaster lineatus Prionospio cirrifera
52.793744 47.863335 46.272692
Polydora ciliata Cumella limicola Monocorophium acherusicum
41.905219 41.327723 39.891356
Oligochaeta Alitta succinea Spio filicornis
36.764590 36.584687 34.878811
Bittium reticulatum Capitella capitata Salvatoria clavata
34.192799 33.198554 32.958829
Loripes orbiculatus Ampelisca diadema Protodorvillea kefersteini
31.347733 30.569510 30.139960
Syllis gracilis Chamelea gallina Heteromastus filiformis
28.542915 28.524575 28.142266
Abra alba Melinna palmata Amphibalanus improvisus
27.461185 26.859546 23.408099
Rissoa splendida Sphaerosyllis hystrix Lagis koreni
17.117581 16.760313 16.471224
Microdeutopus gryllotalpa Leiochone leiopygos Harmothoe reticulata
16.349670 15.331303 15.209809
Tricolia pullus Stenosoma capito Schistomeringos rudolphi
14.667741 13.849706 13.822279
Ampithoe sp. Microphthalmus sp. Tritia reticulata
12.208710 11.854688 11.232600
Polititapes aureus Nemertea Upogebia pusilla
10.998455 10.045834 10.025894
Iphinoe tenella Chironomidae larvae Kellia suborbicularis
9.824385 9.650905 9.077281
Platynereis dumerilii
8.910425
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.2 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.2)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.2))))
)
(plot.top.sp.glms.lvm.zostera.2 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.2,
mapping = aes(x = species, y = log_y_min(count), colour = station),
labs.legend = paste0("Z", as.numeric(unique(abnd.top.sp.glms.lvm.zostera.2$station))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.2 <- lapply(names(glms.lvm.zostera.2.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.2.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera2" etc.
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.2")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.2", "group_"))))
top.sp.abnd.glms.lvm.zostera.2 <- lapply(top.sp.abnd.glms.lvm.zostera.2, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.2 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station == "Poda" ~ 1,
station == "Otmanli" ~ 2,
station == "Vromos" ~ 3,
station == "Gradina" ~ 4,
station == "Ropotamo" ~ 5)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.2.smry <- zoo.abnd.zostera.long.2 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.2 <- map2(top.sp.abnd.glms.lvm.zostera.2, zoo.abnd.zostera.long.2.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.2.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.2 <- map2(top.sp.abnd.glms.lvm.zostera.2, c(8, 8, 4, 8, 4), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
NA
In the case of the seagrasses and case 2 clusters (= stations), the picture is still more unclear.. I suppose this is in no small part because of the differences 2013-14 - very marked for Poda and Otmanli. I suspect the stations changed in these two years (we were looking for Z. noltii in 2014 in particular) - but still, there is much variability. In the future, it’s probably going to be worth it to have more stations in a meadow, if we really want to have an idea of the communities there, and their variability.
The LRs seem to be a bit lower for groups 2, 4, maybe 5 too - still not sure if you can use that as a significance measure.
For now, in group 1 (= Z1), it’s hard to pick some characteristic species - because of the variability between 2013-2014, no doubt. The species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Polydora ciliata, Prionosprio cirrifera (+ others, medium abundance); and the ones with a significantly lower abundance - or even absent - C. gallina, A. ostroumovi, S. clavata, C. limicola, C. costulata, S. hystrix, S. gracilis, T. pullus.
For group 2 (= Z2), the species with higher abundance - which is not really all that high; this group is also loose, hard to distinguish from group 1 - are: S. gracilis, M. lineatus, P. ciliata. The only species with lower abundance - in fact 0 - is Alitta succinea.
For group 3 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A. diadema. The species with lower abundance (or 0) are: B. reticulatum, P. ciliata, P. cirrifera, A. alba, A. succiena, S. clavata, Oligochaeta, A. improvisus.
For group 4 (= Z4), the species with higher abundance are: M. lineatus (very much so); C. limicola, P. kefersteini, C. gallina, C. capitata. The species with lower abundance (or 0) are: P. ciliata, P. cirrifera, A. succinea, A. improvisus, A. alba.
For group 5 (= Z5), the species with higher abundance are: A. ostroumovi, C. capitata, Oligochaeta. The species with lower abundance (or 0) are: R. splendida, T. pullus.
LVM clusters - case 3 Last try: group 1 = Z1-Z2, group 2 = Z3, group 3 = Z4, group 4 = Z5.
Check the model assumptions.
plot(manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3, family = "negative.binomial"))
meanvar.plot(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3, table = TRUE)
START SECTION 2
Plotting if overlay is TRUE
using grouping variable lvm.clusters.zostera.3 160 mean values were 0 and could
not be included in the log-plot
using grouping variable lvm.clusters.zostera.3 160 variance values were 0 and could not
be included in the log-plot
FINISHED SECTION 2
$mean
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 57.6875 50.8125 30.0625 6.25 3.250
2 0.0000 88.0000 0.7500 80.50 2.250
3 238.3750 23.1250 48.6250 24.50 76.375
4 74.2500 124.0000 201.2500 20.25 2.750
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata
1 17.8125 25.375 26.5625
2 2.2500 0.000 0.0000
3 23.0000 0.000 0.2500
4 5.2500 20.000 11.2500
Monocorophium.acherusicum Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi
1 3.9375 7.4375 0.4375 0.0625
2 78.2500 9.5000 0.0000 0.0000
3 1.0000 10.1250 9.8750 3.8750
4 2.2500 4.0000 39.7500 48.2500
Spio.filicornis Microdeutopus.gryllotalpa Abra.alba Cumella.limicola
1 1.4375 2.625 6.3125 0.375
2 38.2500 4.250 0.0000 0.000
3 1.2500 1.875 0.7500 10.250
4 0.0000 14.000 1.7500 6.500
Loripes.orbiculatus Parvicardium.exiguum Protodorvillea.kefersteini
1 1.125 2.8125 0.750
2 0.500 5.5000 0.000
3 8.750 2.8750 9.875
4 4.500 1.5000 1.250
Platynereis.dumerilii Nemertea Syllis.gracilis Alitta.succinea Amphibalanus.improvisus
1 1.5625 2.375 1.625 3.5625 3.3125
2 9.2500 0.500 0.000 0.0000 0.2500
3 1.6250 3.125 0.875 0.0000 0.2500
4 3.7500 1.500 7.250 1.0000 0.5000
Stenosoma.capito Lagis.koreni Schistomeringos.rudolphi Salvatoria.clavata
1 1.500 2.125 1.6875 0.000
2 0.000 0.000 0.0000 1.500
3 0.875 1.125 1.2500 0.375
4 3.750 0.500 0.0000 6.250
Leiochone.leiopygos Melinna.palmata Microphthalmus.sp. Kellia.suborbicularis
1 0.625 0.875 0.000 0.3125
2 0.000 2.750 0.000 0.0000
3 1.875 0.000 0.625 2.3750
4 0.000 0.000 5.000 0.0000
Nototropis.guttatus Chamelea.gallina Perinereis.cultrifera Sphaerosyllis.hystrix
1 0.750 0.0625 0.4375 0.1875
2 0.000 0.0000 0.0000 0.0000
3 0.125 2.0000 0.6250 1.3750
4 2.500 1.2500 2.0000 1.2500
Ampithoe.sp. Harmothoe.reticulata Phoronida Perioculodes.longimanus Rissoa.splendida
1 0.1875 0.750 0.5625 0.375 0.00
2 2.7500 0.000 0.0000 1.250 1.00
3 0.5000 0.125 0.7500 0.125 0.25
4 0.0000 1.000 0.2500 0.750 2.25
Diogenes.pugilator Iphinoe.tenella Platyhelminthes Exogone.naidina
1 0.375 0.6875 0.6875 0.0625
2 0.750 0.0000 0.0000 2.2500
3 0.250 0.2500 0.2500 0.2500
4 0.750 0.0000 0.0000 0.0000
Genetyllis.tuberculata Parthenina.interstincta Tritia.reticulata Cytharella.costulata
1 0.500 0.625 0.5625 0.1875
2 0.250 0.000 0.0000 0.2500
3 0.375 0.250 0.2500 0.6250
4 0.000 0.000 0.2500 0.5000
Syllis.hyalina Tricolia.pullus Turbellaria Harmothoe.imbricata Hydrobia.sp.
1 0.0 0.0625 0.375 0.3125 0.375
2 0.0 0.0000 0.250 0.0000 0.000
3 0.0 0.2500 0.375 0.2500 0.125
4 2.5 1.7500 0.000 0.0000 0.000
Lucinella.divaricata Nereis.pulsatoria Polititapes.aureus Upogebia.pusilla Glycera.sp.
1 0.000 0.4375 0.375 0.00 0.375
2 0.000 0.0000 0.000 0.00 0.000
3 0.875 0.0000 0.125 0.00 0.000
4 0.000 0.0000 0.000 1.75 0.000
Microphthalmus.fragilis Eunice.vittata Glycera.tridactyla Tritia.neritea
1 0.375 0.1875 0.3125 0.250
2 0.000 0.0000 0.0000 0.000
3 0.000 0.0000 0.0000 0.125
4 0.000 0.5000 0.0000 0.000
Chironomidae.larvae Hydrobia.acuta Micromaldane.ornithochaeta Cumella.pygmaea
1 0 0.25 0.000 0.000
2 0 0.00 0.250 0.000
3 0 0.00 0.375 0.375
4 1 0.00 0.000 0.000
Eurydice.dollfusi Hirudinea Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis
1 0.000 0.125 0.0625 0.00 0.0625
2 0.000 0.000 0.0000 0.00 0.0000
3 0.375 0.125 0.1250 0.25 0.2500
4 0.000 0.000 0.2500 0.25 0.0000
Brachynotus.sexdentatus Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa
1 0.0625 0.00 0.0625 0.0625
2 0.0000 0.00 0.0000 0.0000
3 0.1250 0.25 0.0000 0.1250
4 0.0000 0.00 0.2500 0.0000
Abra.sp. Actiniaria Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava
1 0.00 0.0625 0.0625 0.000 0.00
2 0.25 0.0000 0.0000 0.000 0.25
3 0.00 0.0000 0.0000 0.125 0.00
4 0.00 0.0000 0.0000 0.000 0.00
Gastrosaccus.sanctus Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis
1 0.000 0.0625 0.000 0.000
2 0.000 0.0000 0.000 0.000
3 0.125 0.0000 0.125 0.125
4 0.000 0.0000 0.000 0.000
Maldane.glebifex Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.0625 0.0625 0.0625 0.000
2 0.0000 0.0000 0.0000 0.000
3 0.0000 0.0000 0.0000 0.125
4 0.0000 0.0000 0.0000 0.000
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.0625 0.00 0.00 0.00
2 0.0000 0.25 0.25 0.25
3 0.0000 0.00 0.00 0.00
4 0.0000 0.00 0.00 0.00
Thracia.phaseolina
1 0.000
2 0.000
3 0.125
4 0.000
$var
Bittium.reticulatum Capitella.minima Oligochaeta Ampelisca.diadema Mytilaster.lineatus
1 2642.229 2899.2292 700.8625 45.53333 35.000000
2 0.000 3094.0000 2.2500 395.00000 1.583333
3 25180.554 452.6964 3778.2679 735.14286 877.125000
4 3278.250 3103.3333 18570.9167 92.91667 2.916667
Heteromastus.filiformis Prionospio.cirrifera Polydora.ciliata
1 177.495833 919.5833 1000.2625000
2 1.583333 0.0000 0.0000000
3 136.000000 0.0000 0.2142857
4 32.916667 183.3333 34.2500000
Monocorophium.acherusicum Rissoa.membranacea Capitella.capitata Apseudopsis.ostroumovi
1 17.929167 271.729167 1.595833 0.062500
2 3980.916667 11.666667 0.000000 0.000000
3 1.428571 95.553571 82.982143 9.839286
4 4.250000 6.666667 566.250000 196.250000
Spio.filicornis Microdeutopus.gryllotalpa Abra.alba Cumella.limicola
1 2.529167 10.783333 24.362500 1.183333
2 306.250000 16.250000 0.000000 0.000000
3 3.357143 2.696429 1.071429 30.500000
4 0.000000 40.666667 1.583333 23.000000
Loripes.orbiculatus Parvicardium.exiguum Protodorvillea.kefersteini
1 2.383333 18.829167 2.60000
2 1.000000 3.666667 0.00000
3 16.500000 8.410714 16.69643
4 7.000000 1.666667 2.25000
Platynereis.dumerilii Nemertea Syllis.gracilis Alitta.succinea
1 6.929167 12.2500000 7.450000 30.1291667
2 32.916667 0.3333333 0.000000 0.0000000
3 2.553571 3.8392857 2.982143 0.0000000
4 10.916667 3.6666667 22.916667 0.6666667
Amphibalanus.improvisus Stenosoma.capito Lagis.koreni Schistomeringos.rudolphi
1 23.4291667 19.7333333 6.5166667 11.295833
2 0.2500000 0.0000000 0.0000000 0.000000
3 0.2142857 0.9821429 0.9821429 3.642857
4 1.0000000 12.9166667 0.3333333 0.000000
Salvatoria.clavata Leiochone.leiopygos Melinna.palmata Microphthalmus.sp.
1 0.0000000 0.650000 1.183333 0.000
2 0.3333333 0.000000 0.250000 0.000
3 1.1250000 3.839286 0.000000 3.125
4 14.9166667 0.000000 0.000000 22.000
Kellia.suborbicularis Nototropis.guttatus Chamelea.gallina Perinereis.cultrifera
1 0.6291667 1.666667 0.062500 2.2625000
2 0.0000000 0.000000 0.000000 0.0000000
3 13.6964286 0.125000 1.714286 0.8392857
4 0.0000000 14.333333 1.583333 3.3333333
Sphaerosyllis.hystrix Ampithoe.sp. Harmothoe.reticulata Phoronida
1 0.2958333 0.2958333 0.7333333 0.5291667
2 0.0000000 11.5833333 0.0000000 0.0000000
3 2.2678571 0.5714286 0.1250000 1.0714286
4 0.9166667 0.0000000 2.0000000 0.2500000
Perioculodes.longimanus Rissoa.splendida Diogenes.pugilator Iphinoe.tenella
1 0.5166667 0.0000000 0.3833333 2.229167
2 3.5833333 2.0000000 0.2500000 0.000000
3 0.1250000 0.2142857 0.2142857 0.500000
4 0.9166667 6.9166667 0.9166667 0.000000
Platyhelminthes Exogone.naidina Genetyllis.tuberculata Parthenina.interstincta
1 1.4291667 0.0625 2.2666667 2.65
2 0.0000000 8.2500 0.2500000 0.00
3 0.2142857 0.5000 0.5535714 0.50
4 0.0000000 0.0000 0.0000000 0.00
Tritia.reticulata Cytharella.costulata Syllis.hyalina Tricolia.pullus Turbellaria
1 1.4625000 0.1625000 0 0.0625000 1.050
2 0.0000000 0.2500000 0 0.0000000 0.250
3 0.2142857 0.8392857 0 0.2142857 1.125
4 0.2500000 0.3333333 25 2.9166667 0.000
Harmothoe.imbricata Hydrobia.sp. Lucinella.divaricata Nereis.pulsatoria
1 0.6291667 1.183333 0.000000 1.4625
2 0.0000000 0.000000 0.000000 0.0000
3 0.2142857 0.125000 3.267857 0.0000
4 0.0000000 0.000000 0.000000 0.0000
Polititapes.aureus Upogebia.pusilla Glycera.sp. Microphthalmus.fragilis Eunice.vittata
1 0.650 0.000000 0.3833333 2.25 0.2958333
2 0.000 0.000000 0.0000000 0.00 0.0000000
3 0.125 0.000000 0.0000000 0.00 0.0000000
4 0.000 5.583333 0.0000000 0.00 1.0000000
Glycera.tridactyla Tritia.neritea Chironomidae.larvae Hydrobia.acuta
1 0.6291667 0.600 0 1
2 0.0000000 0.000 0 0
3 0.0000000 0.125 0 0
4 0.0000000 0.000 2 0
Micromaldane.ornithochaeta Cumella.pygmaea Eurydice.dollfusi Hirudinea
1 0.0000000 0.0000000 0.0000000 0.1166667
2 0.2500000 0.0000000 0.0000000 0.0000000
3 0.5535714 0.5535714 0.5535714 0.1250000
4 0.0000000 0.0000000 0.0000000 0.0000000
Pseudocuma.longicorne Spisula.subtruncata Tellina.tenuis Brachynotus.sexdentatus
1 0.0625 0.0000000 0.0625000 0.0625
2 0.0000 0.0000000 0.0000000 0.0000
3 0.1250 0.2142857 0.2142857 0.1250
4 0.2500 0.2500000 0.0000000 0.0000
Lentidium.mediterraneum Micronephthys.stammeri Nephtys.cirrosa Abra.sp. Actiniaria
1 0.0000000 0.0625 0.0625 0.00 0.0625
2 0.0000000 0.0000 0.0000 0.25 0.0000
3 0.2142857 0.0000 0.1250 0.00 0.0000
4 0.0000000 0.2500 0.0000 0.00 0.0000
Anadara.kagoshimensis Apherusa.bispinosa Eteone.flava Gastrosaccus.sanctus
1 0.0625 0.000 0.00 0.000
2 0.0000 0.000 0.25 0.000
3 0.0000 0.125 0.00 0.125
4 0.0000 0.000 0.00 0.000
Glycera.unicornis Lepidochitona.cinerea Magelona.papillicornis Maldane.glebifex
1 0.0625 0.000 0.000 0.0625
2 0.0000 0.000 0.000 0.0000
3 0.0000 0.125 0.125 0.0000
4 0.0000 0.000 0.000 0.0000
Mytilus.galloprovincialis Nephtys.kersivalensis Nereis.perivisceralis
1 0.0625 0.0625 0.000
2 0.0000 0.0000 0.000
3 0.0000 0.0000 0.125
4 0.0000 0.0000 0.000
Paradoneis.harpagonea Phyllodoce.sp. Polychaeta.larvae Polygordius.neapolitanus
1 0.0625 0.00 0.00 0.00
2 0.0000 0.25 0.25 0.25
3 0.0000 0.00 0.00 0.00
4 0.0000 0.00 0.00 0.00
Thracia.phaseolina
1 0.000
2 0.000
3 0.125
4 0.000
More or less the same as case 2 before it.
Everything looks more or less fine; fit the model.
glms.lvm.zostera.3 <- manyglm(zoo.mvabnd.zostera ~ lvm.clusters.zostera.3,
family = "negative.binomial")
Explore the fit (residuals, diagnostic plots, etc.).
## residuals vs fitted values
plot(glms.lvm.zostera.3)
## all traditional (g)lm diagnostic plots
plot.manyglm(glms.lvm.zostera.3, which = 1:3)
# ### source mvabund GLM plotting functions modified to use a grey palette - I just can't redo these plots on my own, the function is doing too complicated things internally to scale the x and y axes
# source(here(functions.dir, "default.plot.manyglm_grey.R"))
# source(here(functions.dir, "plot.manyglm_grey.R"))
#
# par(mfrow = c(3,3))
# lapply(3:3, function(i) plot.manyglm.grey(glms.lvm.zostera, which = i, sub.caption = ""))
# par(mfrow = c(3, 3))
Save the model!
write_rds(glms.lvm.zostera.3,
here(save.dir, "glms_lvm_zostera_3.RDS"))
Let’s see the model summary (NB takes a LOT of time if there are many resamplings!).
(glms.lvm.zostera.3.summary <- summary(glms.lvm.zostera.3,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations if you just want to check it out
show.time = "all")
)
Resampling run 0 finished. Time elapsed: 0.00 min ...
Resampling run 100 finished. Time elapsed: 0.36 min ...
Resampling run 200 finished. Time elapsed: 0.72 min ...
Resampling run 300 finished. Time elapsed: 1.07 min ...
Resampling run 400 finished. Time elapsed: 1.43 min ...
Resampling run 500 finished. Time elapsed: 1.78 min ...
Resampling run 600 finished. Time elapsed: 2.14 min ...
Resampling run 700 finished. Time elapsed: 2.50 min ...
Resampling run 800 finished. Time elapsed: 2.86 min ...
Resampling run 900 finished. Time elapsed: 3.21 min ...
Time elapsed: 0 hr 3 min 34 sec
Test statistics:
LR value Pr(>LR)
(Intercept) 1394.5 0.001 ***
lvm.clusters.zostera.32 358.1 0.001 ***
lvm.clusters.zostera.33 423.5 0.001 ***
lvm.clusters.zostera.34 335.3 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate test statistic:
(Intercept) lvm.clusters.zostera.32
LR value Pr(>LR) LR value Pr(>LR)
Abra.alba 47.727 0.001 18.148 0.004
Abra.sp. 6.107 0.547 3.219 0.632
Actiniaria 2.833 0.939 0.446 0.990
Alitta.succinea 8.284 0.276 5.558 0.259
Ampelisca.diadema 45.273 0.001 20.069 0.002
Amphibalanus.improvisus 15.404 0.025 5.997 0.240
Ampithoe.sp. 5.782 0.605 6.490 0.224
Anadara.kagoshimensis 2.833 0.939 0.446 0.990
Apherusa.bispinosa 5.024 0.728 0.000 1.000
Apseudopsis.ostroumovi 21.836 0.019 0.445 0.990
Bittium.reticulatum 159.822 0.001 25.436 0.002
Brachynotus.sexdentatus 4.881 0.753 0.446 0.990
Capitella.capitata 3.390 0.925 2.606 0.762
Capitella.minima 158.885 0.001 0.860 0.984
Chamelea.gallina 21.471 0.019 0.446 0.990
Chironomidae.larvae 9.401 0.220 0.000 0.996
Cumella.limicola 6.366 0.532 2.538 0.774
Cumella.pygmaea 7.559 0.402 0.000 1.000
Cytharella.costulata 10.276 0.177 0.059 0.996
Diogenes.pugilator 7.671 0.396 0.877 0.984
Eteone.flava 6.107 0.547 3.219 0.632
Eunice.vittata 3.020 0.926 1.027 0.968
Eurydice.dollfusi 7.559 0.402 0.000 1.000
Exogone.naidina 6.206 0.534 5.880 0.240
Gastrosaccus.sanctus 5.024 0.728 0.000 1.000
Genetyllis.tuberculata 0.996 0.986 0.185 0.996
Glycera.sp. 4.780 0.753 2.664 0.745
Glycera.tridactyla 2.091 0.973 1.471 0.930
Glycera.unicornis 2.833 0.939 0.446 0.990
Harmothoe.imbricata 3.139 0.926 1.743 0.883
Harmothoe.reticulata 0.997 0.986 5.096 0.313
Heteromastus.filiformis 123.751 0.001 11.936 0.030
Hirudinea 5.835 0.594 0.893 0.984
Hydrobia.acuta 0.471 0.986 0.467 0.989
Hydrobia.sp. 0.941 0.986 1.166 0.950
Iphinoe.tenella 0.542 0.986 2.867 0.704
Kellia.suborbicularis 2.897 0.939 1.598 0.914
Lagis.koreni 6.615 0.510 8.571 0.108
Leiochone.leiopygos 1.958 0.977 4.059 0.486
Lentidium.mediterraneum 9.353 0.224 0.000 0.996
Lepidochitona.cinerea 5.024 0.728 0.000 1.000
Loripes.orbiculatus 0.205 0.986 1.283 0.950
Lucinella.divaricata 5.816 0.600 0.000 1.000
Magelona.papillicornis 5.024 0.728 0.000 1.000
Maldane.glebifex 2.833 0.939 0.446 0.990
Melinna.palmata 0.261 0.986 6.839 0.200
Microdeutopus.gryllotalpa 13.341 0.044 0.870 0.984
Micromaldane.ornithochaeta 9.788 0.198 3.053 0.656
Micronephthys.stammeri 5.710 0.605 0.446 0.990
Microphthalmus.fragilis 0.234 0.986 0.467 0.989
Microphthalmus.sp. 9.420 0.219 0.000 1.000
Monocorophium.acherusicum 26.934 0.005 26.810 0.001
Mytilaster.lineatus 17.441 0.019 0.324 0.993
Mytilus.galloprovincialis 2.833 0.939 0.446 0.990
Nemertea 9.059 0.233 3.320 0.631
Nephtys.cirrosa 4.881 0.753 0.446 0.990
Nephtys.kersivalensis 2.833 0.939 0.446 0.990
Nereis.perivisceralis 5.024 0.728 0.000 1.000
Nereis.pulsatoria 0.491 0.986 0.980 0.971
Nototropis.guttatus 0.323 0.986 2.976 0.669
Oligochaeta 126.966 0.001 15.119 0.007
Paradoneis.harpagonea 2.833 0.939 0.446 0.990
Parthenina.interstincta 0.301 0.986 1.536 0.924
Parvicardium.exiguum 14.018 0.039 1.502 0.930
Perinereis.cultrifera 1.844 0.981 2.060 0.867
Perioculodes.longimanus 3.981 0.870 2.049 0.869
Phoronida 3.161 0.926 3.937 0.509
Phyllodoce.sp. 6.107 0.547 3.219 0.632
Platyhelminthes 0.644 0.986 3.118 0.642
Platynereis.dumerilii 1.824 0.981 7.414 0.173
Polititapes.aureus 2.789 0.939 2.081 0.863
Polychaeta.larvae 6.107 0.547 3.219 0.632
Polydora.ciliata 81.609 0.001 15.796 0.005
Polygordius.neapolitanus 6.107 0.547 3.219 0.632
Prionospio.cirrifera 91.317 0.001 19.403 0.002
Protodorvillea.kefersteini 0.820 0.986 4.746 0.362
Pseudocuma.longicorne 7.101 0.453 0.446 0.990
Rissoa.membranacea 50.503 0.001 0.160 0.996
Rissoa.splendida 17.291 0.019 10.624 0.047
Salvatoria.clavata 25.703 0.006 18.216 0.004
Schistomeringos.rudolphi 0.703 0.986 2.760 0.729
Sphaerosyllis.hystrix 10.217 0.179 1.307 0.949
Spio.filicornis 1.763 0.981 24.367 0.002
Spisula.subtruncata 11.731 0.089 0.000 1.000
Stenosoma.capito 1.007 0.986 4.916 0.344
Syllis.gracilis 1.598 0.981 5.347 0.283
Syllis.hyalina 3.627 0.894 0.000 0.996
Tellina.tenuis 7.795 0.377 0.446 0.990
Thracia.phaseolina 5.024 0.728 0.000 1.000
Tricolia.pullus 13.409 0.043 0.445 0.990
Tritia.neritea 1.950 0.978 1.104 0.961
Tritia.reticulata 1.128 0.983 2.531 0.777
Turbellaria 1.449 0.981 0.056 0.996
Upogebia.pusilla 8.660 0.248 0.000 1.000
lvm.clusters.zostera.33 lvm.clusters.zostera.34
LR value Pr(>LR) LR value
Abra.alba 15.271 0.010 4.994
Abra.sp. 0.000 1.000 0.000
Actiniaria 0.811 1.000 0.446
Alitta.succinea 9.839 0.093 1.072
Ampelisca.diadema 8.950 0.133 4.744
Amphibalanus.improvisus 10.321 0.082 4.259
Ampithoe.sp. 1.064 1.000 1.175
Anadara.kagoshimensis 0.811 1.000 0.446
Apherusa.bispinosa 2.197 0.995 0.000
Apseudopsis.ostroumovi 32.536 0.001 46.727
Bittium.reticulatum 10.139 0.084 0.217
Brachynotus.sexdentatus 0.236 1.000 0.446
Capitella.capitata 18.553 0.002 25.101
Capitella.minima 2.320 0.991 2.355
Chamelea.gallina 22.820 0.002 11.133
Chironomidae.larvae 0.000 1.000 7.676
Cumella.limicola 27.637 0.001 20.452
Cumella.pygmaea 4.909 0.612 0.000
Cytharella.costulata 2.833 0.963 1.046
Diogenes.pugilator 0.263 1.000 0.877
Eteone.flava 0.000 1.000 0.000
Eunice.vittata 1.873 0.999 0.483
Eurydice.dollfusi 4.909 0.612 0.000
Exogone.naidina 1.012 1.000 0.416
Gastrosaccus.sanctus 2.197 0.995 0.000
Genetyllis.tuberculata 0.063 1.000 1.888
Glycera.sp. 4.773 0.632 2.664
Glycera.tridactyla 2.652 0.977 1.471
Glycera.unicornis 0.811 1.000 0.446
Harmothoe.imbricata 0.047 1.000 1.743
Harmothoe.reticulata 4.553 0.677 0.222
Heteromastus.filiformis 0.653 1.000 5.770
Hirudinea 0.000 1.000 0.893
Hydrobia.acuta 0.845 1.000 0.467
Hydrobia.sp. 0.396 1.000 1.166
Iphinoe.tenella 0.979 1.000 2.867
Kellia.suborbicularis 4.332 0.730 1.598
Lagis.koreni 1.378 1.000 2.889
Leiochone.leiopygos 4.772 0.633 4.059
Lentidium.mediterraneum 4.394 0.718 0.000
Lepidochitona.cinerea 2.197 0.995 0.000
Loripes.orbiculatus 23.881 0.002 9.848
Lucinella.divaricata 4.958 0.603 0.000
Magelona.papillicornis 2.197 0.995 0.000
Maldane.glebifex 0.811 1.000 0.446
Melinna.palmata 11.353 0.053 6.248
Microdeutopus.gryllotalpa 0.549 1.000 10.829
Micromaldane.ornithochaeta 5.601 0.498 0.000
Micronephthys.stammeri 0.811 1.000 0.893
Microphthalmus.fragilis 0.845 1.000 0.467
Microphthalmus.sp. 6.434 0.390 10.604
Monocorophium.acherusicum 6.554 0.382 0.880
Mytilaster.lineatus 28.164 0.001 0.072
Mytilus.galloprovincialis 0.811 1.000 0.446
Nemertea 0.344 1.000 0.460
Nephtys.cirrosa 0.236 1.000 0.446
Nephtys.kersivalensis 0.811 1.000 0.446
Nereis.perivisceralis 2.197 0.995 0.000
Nereis.pulsatoria 1.766 0.999 0.980
Nototropis.guttatus 2.256 0.995 1.529
Oligochaeta 1.311 1.000 12.946
Paradoneis.harpagonea 0.811 1.000 0.446
Parthenina.interstincta 0.348 1.000 1.536
Parvicardium.exiguum 0.002 1.000 0.885
Perinereis.cultrifera 0.146 1.000 2.073
Perioculodes.longimanus 1.060 1.000 0.584
Phoronida 0.286 1.000 0.722
Phyllodoce.sp. 0.000 1.000 0.000
Platyhelminthes 1.116 1.000 3.118
Platynereis.dumerilii 0.005 1.000 1.755
Polititapes.aureus 0.920 1.000 2.081
Polychaeta.larvae 0.000 1.000 0.000
Polydora.ciliata 20.565 0.002 1.275
Polygordius.neapolitanus 0.000 1.000 0.000
Prionospio.cirrifera 31.954 0.001 0.160
Protodorvillea.kefersteini 18.752 0.002 0.636
Pseudocuma.longicorne 0.236 1.000 0.893
Rissoa.membranacea 0.416 1.000 0.818
Rissoa.splendida 4.195 0.748 16.268
Salvatoria.clavata 6.518 0.382 31.870
Schistomeringos.rudolphi 0.066 1.000 2.760
Sphaerosyllis.hystrix 8.440 0.170 5.462
Spio.filicornis 0.088 1.000 7.758
Spisula.subtruncata 4.394 0.719 3.219
Stenosoma.capito 0.492 1.000 1.211
Syllis.gracilis 0.703 1.000 3.662
Syllis.hyalina 0.000 1.000 3.584
Tellina.tenuis 1.386 1.000 0.446
Thracia.phaseolina 2.197 0.995 0.000
Tricolia.pullus 1.365 1.000 11.315
Tritia.neritea 0.205 1.000 1.104
Tritia.reticulata 0.624 1.000 0.361
Turbellaria 0.000 1.000 1.475
Upogebia.pusilla 0.000 1.000 8.030
Pr(>LR)
Abra.alba 0.596
Abra.sp. 1.000
Actiniaria 1.000
Alitta.succinea 1.000
Ampelisca.diadema 0.632
Amphibalanus.improvisus 0.702
Ampithoe.sp. 1.000
Anadara.kagoshimensis 1.000
Apherusa.bispinosa 1.000
Apseudopsis.ostroumovi 0.001
Bittium.reticulatum 1.000
Brachynotus.sexdentatus 1.000
Capitella.capitata 0.001
Capitella.minima 0.965
Chamelea.gallina 0.068
Chironomidae.larvae 0.237
Cumella.limicola 0.001
Cumella.pygmaea 1.000
Cytharella.costulata 1.000
Diogenes.pugilator 1.000
Eteone.flava 1.000
Eunice.vittata 1.000
Eurydice.dollfusi 1.000
Exogone.naidina 1.000
Gastrosaccus.sanctus 1.000
Genetyllis.tuberculata 0.988
Glycera.sp. 0.944
Glycera.tridactyla 0.997
Glycera.unicornis 1.000
Harmothoe.imbricata 0.991
Harmothoe.reticulata 1.000
Heteromastus.filiformis 0.459
Hirudinea 1.000
Hydrobia.acuta 1.000
Hydrobia.sp. 1.000
Iphinoe.tenella 0.903
Kellia.suborbicularis 0.995
Lagis.koreni 0.903
Leiochone.leiopygos 0.735
Lentidium.mediterraneum 1.000
Lepidochitona.cinerea 1.000
Loripes.orbiculatus 0.104
Lucinella.divaricata 1.000
Magelona.papillicornis 1.000
Maldane.glebifex 1.000
Melinna.palmata 0.399
Microdeutopus.gryllotalpa 0.075
Micromaldane.ornithochaeta 1.000
Micronephthys.stammeri 1.000
Microphthalmus.fragilis 1.000
Microphthalmus.sp. 0.079
Monocorophium.acherusicum 1.000
Mytilaster.lineatus 1.000
Mytilus.galloprovincialis 1.000
Nemertea 1.000
Nephtys.cirrosa 1.000
Nephtys.kersivalensis 1.000
Nereis.perivisceralis 1.000
Nereis.pulsatoria 1.000
Nototropis.guttatus 0.996
Oligochaeta 0.027
Paradoneis.harpagonea 1.000
Parthenina.interstincta 0.996
Parvicardium.exiguum 1.000
Perinereis.cultrifera 0.980
Perioculodes.longimanus 1.000
Phoronida 1.000
Phyllodoce.sp. 1.000
Platyhelminthes 0.867
Platynereis.dumerilii 0.991
Polititapes.aureus 0.979
Polychaeta.larvae 1.000
Polydora.ciliata 1.000
Polygordius.neapolitanus 1.000
Prionospio.cirrifera 1.000
Protodorvillea.kefersteini 1.000
Pseudocuma.longicorne 1.000
Rissoa.membranacea 1.000
Rissoa.splendida 0.005
Salvatoria.clavata 0.001
Schistomeringos.rudolphi 0.926
Sphaerosyllis.hystrix 0.498
Spio.filicornis 0.233
Spisula.subtruncata 0.854
Stenosoma.capito 1.000
Syllis.gracilis 0.797
Syllis.hyalina 0.797
Tellina.tenuis 1.000
Thracia.phaseolina 1.000
Tricolia.pullus 0.060
Tritia.neritea 1.000
Tritia.reticulata 1.000
Turbellaria 0.997
Upogebia.pusilla 0.216
Arguments: with 999 resampling iterations using pit.trap resampling and response assumed to be uncorrelated
Likelihood Ratio statistic: 965.7, p-value: 0.001
Univariate test statistic:
Abra.alba Abra.sp. Actiniaria Alitta.succinea Ampelisca.diadema
LR value 27.438 4.159 1.386 13.528 20.153
Pr(>LR) 0.002 0.933 0.994 0.119 0.010
Amphibalanus.improvisus Ampithoe.sp. Anadara.kagoshimensis Apherusa.bispinosa
LR value 15.231 8.727 1.386 2.773
Pr(>LR) 0.076 0.517 0.994 0.983
Apseudopsis.ostroumovi Bittium.reticulatum Brachynotus.sexdentatus
LR value 51.414 34.189 1.386
Pr(>LR) 0.001 0.001 0.994
Capitella.capitata Capitella.minima Chamelea.gallina Chironomidae.larvae
LR value 30.336 6.793 27.138 9.651
Pr(>LR) 0.002 0.718 0.002 0.396
Cumella.limicola Cumella.pygmaea Cytharella.costulata Diogenes.pugilator
LR value 33.634 6.118 3.195 2.380
Pr(>LR) 0.001 0.804 0.983 0.991
Eteone.flava Eunice.vittata Eurydice.dollfusi Exogone.naidina
LR value 4.159 3.572 6.118 7.227
Pr(>LR) 0.933 0.977 0.804 0.690
Gastrosaccus.sanctus Genetyllis.tuberculata Glycera.sp. Glycera.tridactyla
LR value 2.773 1.966 7.992 4.484
Pr(>LR) 0.983 0.994 0.575 0.925
Glycera.unicornis Harmothoe.imbricata Harmothoe.reticulata
LR value 1.386 3.192 9.416
Pr(>LR) 0.994 0.983 0.421
Heteromastus.filiformis Hirudinea Hydrobia.acuta Hydrobia.sp. Iphinoe.tenella
LR value 17.463 1.726 1.438 2.250 5.470
Pr(>LR) 0.027 0.994 0.994 0.991 0.869
Kellia.suborbicularis Lagis.koreni Leiochone.leiopygos Lentidium.mediterraneum
LR value 8.973 10.472 14.995 5.545
Pr(>LR) 0.469 0.312 0.077 0.865
Lepidochitona.cinerea Loripes.orbiculatus Lucinella.divaricata
LR value 2.773 28.176 6.172
Pr(>LR) 0.983 0.002 0.795
Magelona.papillicornis Maldane.glebifex Melinna.palmata
LR value 2.773 1.386 24.203
Pr(>LR) 0.983 0.994 0.003
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta Micronephthys.stammeri
LR value 13.747 6.655 2.773
Pr(>LR) 0.113 0.721 0.983
Microphthalmus.fragilis Microphthalmus.sp. Monocorophium.acherusicum
LR value 1.437 11.855 39.678
Pr(>LR) 0.994 0.208 0.001
Mytilaster.lineatus Mytilus.galloprovincialis Nemertea Nephtys.cirrosa
LR value 34.910 1.386 4.566 1.386
Pr(>LR) 0.001 0.994 0.925 0.994
Nephtys.kersivalensis Nereis.perivisceralis Nereis.pulsatoria
LR value 1.386 2.773 2.987
Pr(>LR) 0.994 0.983 0.983
Nototropis.guttatus Oligochaeta Paradoneis.harpagonea Parthenina.interstincta
LR value 7.095 27.378 1.386 2.884
Pr(>LR) 0.690 0.002 0.994 0.983
Parvicardium.exiguum Perinereis.cultrifera Perioculodes.longimanus Phoronida
LR value 2.911 4.654 4.365 5.372
Pr(>LR) 0.983 0.925 0.925 0.879
Phyllodoce.sp. Platyhelminthes Platynereis.dumerilii Polititapes.aureus
LR value 4.159 5.906 8.537 4.065
Pr(>LR) 0.933 0.834 0.523 0.943
Polychaeta.larvae Polydora.ciliata Polygordius.neapolitanus
LR value 4.159 29.489 4.159
Pr(>LR) 0.933 0.002 0.933
Prionospio.cirrifera Protodorvillea.kefersteini Pseudocuma.longicorne
LR value 41.351 25.261 1.726
Pr(>LR) 0.001 0.002 0.994
Rissoa.membranacea Rissoa.splendida Salvatoria.clavata Schistomeringos.rudolphi
LR value 1.701 17.118 32.959 5.024
Pr(>LR) 0.994 0.032 0.001 0.900
Sphaerosyllis.hystrix Spio.filicornis Spisula.subtruncata Stenosoma.capito
LR value 12.736 34.632 5.885 7.189
Pr(>LR) 0.154 0.001 0.834 0.690
Syllis.gracilis Syllis.hyalina Tellina.tenuis Thracia.phaseolina
LR value 10.712 4.555 3.112 2.773
Pr(>LR) 0.306 0.925 0.983 0.983
Tricolia.pullus Tritia.neritea Tritia.reticulata Turbellaria Upogebia.pusilla
LR value 13.288 2.037 2.932 1.527 10.026
Pr(>LR) 0.128 0.994 0.983 0.994 0.353
Arguments:
Test statistics calculated assuming response assumed to be uncorrelated
P-value calculated using 999 resampling iterations via pit.trap resampling (to account for correlation in testing).
The factor is highly significant according to the models.
Again, save the summary for safekeeping, but also run an anova.
write_rds(glms.lvm.zostera.3.summary,
here(save.dir, "glms_lvm_zostera_3_summary.RDS"))
Run the anova on the model.
(glms.lvm.zostera.3.aov <- anova.manyglm(glms.lvm.zostera.3,
test = "LR", p.uni = "adjusted",
nBoot = 999, ## limit the number of permutations for a shorter run time
show.time = "all")
)
Resampling begins for test 1.
Resampling run 0 finished. Time elapsed: 0.00 minutes...
Resampling run 100 finished. Time elapsed: 0.12 minutes...
Resampling run 200 finished. Time elapsed: 0.24 minutes...
Resampling run 300 finished. Time elapsed: 0.36 minutes...
Resampling run 400 finished. Time elapsed: 0.48 minutes...
Resampling run 500 finished. Time elapsed: 0.59 minutes...
Resampling run 600 finished. Time elapsed: 0.71 minutes...
Resampling run 700 finished. Time elapsed: 0.83 minutes...
Resampling run 800 finished. Time elapsed: 0.95 minutes...
Resampling run 900 finished. Time elapsed: 1.08 minutes...
Time elapsed: 0 hr 1 min 11 sec
Analysis of Deviance Table
Model: manyglm(formula = zoo.mvabnd.zostera ~ lvm.clusters.zostera.3,
Model: family = "negative.binomial")
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
(Intercept) 31
lvm.clusters.zostera.3 28 3 965.7 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 27.438 0.004 4.159 1.000 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 13.528 0.176 20.153 0.018
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 15.231 0.106 8.727 0.672
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 2.773 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 51.414 0.001 34.189 0.002
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 30.336 0.002
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 6.793 0.895 27.138 0.004
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 9.651 0.562 33.634 0.002
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 6.118 0.944 3.195 1.000
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 2.38 1.000 4.159 1.000 3.572
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 6.118 0.944 7.227 0.875
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 1.966 1.000
Glycera.sp. Glycera.tridactyla
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 7.992 0.788 4.484 0.993
Glycera.unicornis Harmothoe.imbricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 3.192 1.000
Harmothoe.reticulata Heteromastus.filiformis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 9.416 0.602 17.463 0.050
Hirudinea Hydrobia.acuta Hydrobia.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.726 1.000 1.438 1.000 2.25 1.000
Iphinoe.tenella Kellia.suborbicularis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.47 0.983 8.973 0.652
Lagis.koreni Leiochone.leiopygos
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 10.472 0.437 14.995 0.108
Lentidium.mediterraneum Lepidochitona.cinerea
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.545 0.981 2.773 1.000
Loripes.orbiculatus Lucinella.divaricata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 28.176 0.004 6.172 0.942
Magelona.papillicornis Maldane.glebifex
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 1.386 1.000
Melinna.palmata Microdeutopus.gryllotalpa
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 24.203 0.004 13.747 0.172
Micromaldane.ornithochaeta Micronephthys.stammeri
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 6.655 0.908 2.773
Microphthalmus.fragilis Microphthalmus.sp.
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 1.000 1.437 1.000 11.855
Monocorophium.acherusicum Mytilaster.lineatus
Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 0.278 39.678 0.001 34.91
Mytilus.galloprovincialis Nemertea
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 0.002 1.386 1.000 4.566 0.991
Nephtys.cirrosa Nephtys.kersivalensis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 1.386 1.000
Nereis.perivisceralis Nereis.pulsatoria
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.773 1.000 2.987 1.000
Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 7.095 0.879 27.378 0.004
Paradoneis.harpagonea Parthenina.interstincta
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.386 1.000 2.884 1.000
Parvicardium.exiguum Perinereis.cultrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 2.911 1.000 4.654 0.991
Perioculodes.longimanus Phoronida
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.365 0.994 5.372 0.983
Phyllodoce.sp. Platyhelminthes
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 5.906 0.957
Platynereis.dumerilii Polititapes.aureus
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 8.537 0.700 4.065 1.000
Polychaeta.larvae Polydora.ciliata
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 29.489 0.002
Polygordius.neapolitanus Prionospio.cirrifera
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 4.159 1.000 41.351 0.001
Protodorvillea.kefersteini Pseudocuma.longicorne
Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 25.261 0.004 1.726
Rissoa.membranacea Rissoa.splendida
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 1.701 1.000 17.118 0.058
Salvatoria.clavata Schistomeringos.rudolphi
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 32.959 0.002 5.024 0.985
Sphaerosyllis.hystrix Spio.filicornis
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 12.736 0.216 34.632 0.002
Spisula.subtruncata Stenosoma.capito
Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 5.885 0.957 7.189 0.879
Syllis.gracilis Syllis.hyalina Tellina.tenuis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 10.712 0.407 4.555 0.992 3.112
Thracia.phaseolina Tricolia.pullus
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 2.773 1.000 13.288 0.176
Tritia.neritea Tritia.reticulata Turbellaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev
(Intercept)
lvm.clusters.zostera.3 2.037 1.000 2.932 1.000 1.527
Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev)
(Intercept)
lvm.clusters.zostera.3 1.000 10.026 0.493
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Save the ANOVA, too.
write_rds(glms.lvm.zostera.3.aov,
here(save.dir, "glms_lvm_zostera_3_anova.RDS"))
NOW let’s get the taxa with the highest contributions to the tested pattern.
## get the top contributing species for the initial zostera GLMs
(top.sp.glms.lvm.zostera.3 <- top_n_sp_glm(glms.lvm.zostera.3.aov, tot.dev.expl = 0.75)
)
[1] "Total deviance explained: 0.82"
Apseudopsis.ostroumovi Prionospio.cirrifera Monocorophium.acherusicum
51.413970 41.350974 39.678240
Mytilaster.lineatus Spio.filicornis Bittium.reticulatum
34.909870 34.632433 34.188511
Cumella.limicola Salvatoria.clavata Capitella.capitata
33.633823 32.958829 30.336457
Polydora.ciliata Loripes.orbiculatus Abra.alba
29.489490 28.175552 27.437828
Oligochaeta Chamelea.gallina Protodorvillea.kefersteini
27.378351 27.138297 25.260955
Melinna.palmata Ampelisca.diadema Heteromastus.filiformis
24.202974 20.153373 17.463258
Rissoa.splendida Amphibalanus.improvisus Leiochone.leiopygos
17.117581 15.231072 14.994968
Microdeutopus.gryllotalpa Alitta.succinea Tricolia.pullus
13.746930 13.527612 13.287785
Sphaerosyllis.hystrix Microphthalmus.sp. Syllis.gracilis
12.736427 11.854688 10.711746
Lagis.koreni Upogebia.pusilla Chironomidae.larvae
10.471907 10.025894 9.650905
Harmothoe.reticulata Kellia.suborbicularis Ampithoe.sp.
9.416160 8.973262 8.726571
Platynereis.dumerilii Glycera.sp. Exogone.naidina
8.536643 7.992056 7.227403
Stenosoma.capito Nototropis.guttatus Capitella.minima
7.189015 7.094867 6.793482
Micromaldane.ornithochaeta
6.655076
## unfortunately, mvabund likes to rename my species when converting the data to matrix (no spaces in names), and since I'm going to look them up in my initial untransformed count data, I have to change them back..
names(top.sp.glms.lvm.zostera.3) <- names(top.sp.glms.lvm.zostera.3) %>%
str_replace(pattern = "\\.", replacement = " ")
top.sp.glms.lvm.zostera.3
Apseudopsis ostroumovi Prionospio cirrifera Monocorophium acherusicum
51.413970 41.350974 39.678240
Mytilaster lineatus Spio filicornis Bittium reticulatum
34.909870 34.632433 34.188511
Cumella limicola Salvatoria clavata Capitella capitata
33.633823 32.958829 30.336457
Polydora ciliata Loripes orbiculatus Abra alba
29.489490 28.175552 27.437828
Oligochaeta Chamelea gallina Protodorvillea kefersteini
27.378351 27.138297 25.260955
Melinna palmata Ampelisca diadema Heteromastus filiformis
24.202974 20.153373 17.463258
Rissoa splendida Amphibalanus improvisus Leiochone leiopygos
17.117581 15.231072 14.994968
Microdeutopus gryllotalpa Alitta succinea Tricolia pullus
13.746930 13.527612 13.287785
Sphaerosyllis hystrix Microphthalmus sp. Syllis gracilis
12.736427 11.854688 10.711746
Lagis koreni Upogebia pusilla Chironomidae larvae
10.471907 10.025894 9.650905
Harmothoe reticulata Kellia suborbicularis Ampithoe sp.
9.416160 8.973262 8.726571
Platynereis dumerilii Glycera sp. Exogone naidina
8.536643 7.992056 7.227403
Stenosoma capito Nototropis guttatus Capitella minima
7.189015 7.094867 6.793482
Micromaldane ornithochaeta
6.655076
Try to plot these top contributing species - for whatever that’s worth, because 50 species on a plot is still a monstrosity.
## get the species and their abundances from the original count data, and transform them to long format
(abnd.top.sp.glms.lvm.zostera.3 <- zoo.abnd.zostera %>%
select(station, names(top.sp.glms.lvm.zostera.3)) %>%
gather(key = "species", value = "count", -station) %>%
## turn species into a factor, or you'll be very very sorry later, when they're out of order on the plot. NB need to be in REVERSE order, because ggplot plots from bottom to top, and I want the top-contributing species on top.
mutate(species = factor(species, levels = rev(names(top.sp.glms.lvm.zostera.3))))
)
(plot.top.sp.glms.lvm.zostera.3 <- plot_top_n(abnd.top.sp.glms.lvm.zostera.3,
mapping = aes(x = species, y = log_y_min(count), colour = station),
labs.legend = paste0("Z", as.numeric(unique(abnd.top.sp.glms.lvm.zostera.3$station))),
lab.y = "Abundance (log(y/min + 1))",
palette = "Set2"
) +
theme(legend.position = "top")
)
Extract the top-contributing species to each cluster (this same nightmare above, but as a table). This chunk is STILL hopelessly ugly and clumsy.
top.sp.abnd.glms.lvm.zostera.3 <- lapply(names(glms.lvm.zostera.3.summary$aliased), function(x) top_sp_glms_table(glms.lvm.zostera.3.summary, x, p = 0.05))
## fix species names (remove dot)
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% mutate(species = str_replace(species, pattern = "\\.", replacement = " ")))
## rename columns (= group names) - right now they are something like "lvm.clusters.zostera3" etc.
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% rename_at(vars(contains("lvm.clusters.zostera.3")), list(~str_replace_all(., pattern = "lvm.clusters.zostera.3", "group_"))))
top.sp.abnd.glms.lvm.zostera.3 <- lapply(top.sp.abnd.glms.lvm.zostera.3, function(x) x %>% rename_at(vars(contains("Intercept")), list(~str_replace_all(., pattern = "\\(Intercept\\)", "group_1"))))
## pull the abundances from the original count df and add to the summary glm tables
## make a long df of abundances & add clusters
zoo.abnd.zostera.long.3 <- zoo.abnd.zostera %>%
select(-c(month:replicate)) %>%
gather(key = "species", value = "count", -station) %>%
mutate(group = case_when(station %in% c("Poda", "Otmanli") ~ 1,
station == "Vromos" ~ 2,
station == "Gradina" ~ 3,
station == "Ropotamo" ~ 4)
)
## sum sp abundances by group; nest by group
zoo.abnd.zostera.long.3.smry <- zoo.abnd.zostera.long.3 %>%
group_by(species, group) %>%
summarise(total_count = sum(count)) %>%
group_by(group) %>%
nest()
## add the counts to the group dfs - wow that's an ugly, ugly hack. Wish I had more time to write this up properly..
top.sp.abnd.glms.lvm.zostera.3 <- map2(top.sp.abnd.glms.lvm.zostera.3, zoo.abnd.zostera.long.3.smry %>% pull(group), ~left_join(.x, zoo.abnd.zostera.long.3.smry %>% filter(group == .y) %>% unnest(), by = "species"))
## since these are sum counts over all the replicates (that's why the monstrous numbers), average them to be mean counts per group. NB different groups consist of different numbers of replicates, b.c. some groups consist of more than one station
(top.sp.abnd.glms.lvm.zostera.3 <- map2(top.sp.abnd.glms.lvm.zostera.3, c(16, 4, 8, 4), function(x, y) x %>% mutate(mean_count = total_count/y))
)
[[1]]
[[2]]
[[3]]
[[4]]
NA
In the case of the seagrasses and case 3 clusters, the picture is still more confusing..
The LRs seem to be a bit lower for groups 2 and 3, maybe 4 too - still not sure if you can use that as a significance measure.
For now, in group 1 (= Z1-Z2), the species/taxa with significantly higher abundance are: Bittium reticulatum, Capitella minima, Oligochaeta, Heteromastus filiformis, Polydora ciliata, Prionosprio cirrifera, Rissoa membranacea, Abra alba, Ampelisca diadema (+ others, medium abundance); and the ones with a significantly lower abundance - or even absent - S. clavata, A. ostroumovi, C. gallina, T. pullus.
There are more species singled out for this cluster, probably because of the variability between the two years of sampling.
For group 2 (= Z3), the species with higher abundance are: M. acherusicum, S. filicornis, A. diadema. The species with lower abundance - in fact 0 - are B. reticulatum, P. cirrifera, S. clavata, A. alba, P. ciliata, oligochaetes, H. filiformis, R. splendida.
For group 3 (= Z4), the species with higher abundance are: M. lineatus, less so - C. limicola, L. orbiculatus, P. kefersteini, C. capitata, etc. The species with lower abundance (or 0) are: P. cirrifera, P. ciliata, A. alba, etc.
For group 4 (= Z5), the species with higher abundance are: A. ostoumovi, C. capitata, oligochaetes (very abundant, but with a small LR - nice!). The species with lower abundance (or 0) are: R. splendida, S. clavata, C. limicola.
All in all, I think that group 1 (Poda + Otmanli) holds, and so does group 2 (Vromos). The question is whether it makes sense to separate Gradina and Ropotamo into 2 groups, or if they make more sense together. Ropotamo is characterized by a very high number of oligochaetes, while Gradina’s most distinguishing characteristic is the high number of M. lineatus - mostly very small ones, attached to the rhizomes, close to the sediment surface I presume. Both stations have C. limicola in medium abundance, which is not present anywhere else.
Try to compare the three models..
glms.lvm.zostera.comp
Analysis of Deviance Table
glms.lvm.zostera.1: zoo.mvabnd.zostera ~ lvm.clusters.zostera.1
glms.lvm.zostera.3: zoo.mvabnd.zostera ~ lvm.clusters.zostera.3
glms.lvm.zostera.2: zoo.mvabnd.zostera ~ lvm.clusters.zostera.2
Multivariate test:
Res.Df Df.diff Dev Pr(>Dev)
glms.lvm.zostera.1 29
glms.lvm.zostera.3 28 1 262.1 0.002 **
glms.lvm.zostera.2 27 1 256.8 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Univariate Tests:
Abra.alba Abra.sp. Actiniaria
Dev Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.506 0.999 0 1.000 0 1.000
glms.lvm.zostera.2 0.023 1.000 0 1.000 1.386 1.000
Alitta.succinea Ampelisca.diadema
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 4.881 0.681 0.094 1.000
glms.lvm.zostera.2 23.057 0.001 10.416 0.115
Amphibalanus.improvisus Ampithoe.sp.
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.349 1.000 2.405 0.970
glms.lvm.zostera.2 8.177 0.274 3.482 0.937
Anadara.kagoshimensis Apherusa.bispinosa
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0.811 1.000
glms.lvm.zostera.2 1.386 1.000 0 1.000
Apseudopsis.ostroumovi Bittium.reticulatum
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 20.559 0.003 3.006 0.943
glms.lvm.zostera.2 1.38 1.000 0.004 1.000
Brachynotus.sexdentatus Capitella.capitata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 4.124 0.797
glms.lvm.zostera.2 1.386 1.000 2.862 0.965
Capitella.minima Chamelea.gallina
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 6.055 0.509 0.908 1.000
glms.lvm.zostera.2 0.433 1.000 1.386 0.996
Chironomidae.larvae Cumella.limicola
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 5.456 0.616 1.089 1.000
glms.lvm.zostera.2 0 1.000 7.694 0.334
Cumella.pygmaea Cytharella.costulata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.891 0.995 0.073 1.000
glms.lvm.zostera.2 0 1.000 4.159 0.859
Diogenes.pugilator Eteone.flava Eunice.vittata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.483 1.000 0 1.000 2.7
glms.lvm.zostera.2 0.68 1.000 0 1.000 0.208
Eurydice.dollfusi Exogone.naidina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.963 1.891 0.995 1.284 1.000
glms.lvm.zostera.2 1.000 0 1.000 1.25 1.000
Gastrosaccus.sanctus Genetyllis.tuberculata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 1.499 0.999
glms.lvm.zostera.2 0 1.000 5.746 0.594
Glycera.sp. Glycera.tridactyla Glycera.unicornis
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0 1.000 0
glms.lvm.zostera.2 0 1.000 1.033 1.000 1.386
Harmothoe.imbricata Harmothoe.reticulata
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.352 1.000 4.286 0.775
glms.lvm.zostera.2 1.000 5.561 0.600 5.794 0.583
Heteromastus.filiformis Hirudinea Hydrobia.acuta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 7.242 0.344 0.811 1.000 0
glms.lvm.zostera.2 10.679 0.104 0 1.000 1.494
Hydrobia.sp. Iphinoe.tenella
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 0.569 1.000 1.31 1.000
glms.lvm.zostera.2 0.996 3.709 0.931 4.354 0.838
Kellia.suborbicularis Lagis.koreni
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 4.389 0.773 0.767 1.000
glms.lvm.zostera.2 0.104 1.000 5.999 0.563
Leiochone.leiopygos Lentidium.mediterraneum
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 8.777 0.191 1.622 0.999
glms.lvm.zostera.2 0.336 1.000 0 1.000
Lepidochitona.cinerea Loripes.orbiculatus
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.811 1.000 2.917 0.946
glms.lvm.zostera.2 0 1.000 3.172 0.942
Lucinella.divaricata Magelona.papillicornis
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.915 0.995 0.811 1.000
glms.lvm.zostera.2 0 1.000 0 1.000
Maldane.glebifex Melinna.palmata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0 1.000
glms.lvm.zostera.2 1.386 0.998 2.657 0.985
Microdeutopus.gryllotalpa Micromaldane.ornithochaeta
Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 10.659 0.084 2.111
glms.lvm.zostera.2 2.603 0.985 0
Micronephthys.stammeri Microphthalmus.fragilis
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0.995 2.197 0.994 0
glms.lvm.zostera.2 1.000 1.386 1.000 1.494
Microphthalmus.sp. Monocorophium.acherusicum
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 2.459 0.970 1.35
glms.lvm.zostera.2 0.996 0 1.000 0.213
Mytilaster.lineatus Mytilus.galloprovincialis
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 12.977 0.028 0
glms.lvm.zostera.2 1.000 12.953 0.043 1.386
Nemertea Nephtys.cirrosa
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.013 1.000 0.811 1.000
glms.lvm.zostera.2 1.000 5.479 0.640 1.386 1.000
Nephtys.kersivalensis Nereis.perivisceralis
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 0.811 1.000
glms.lvm.zostera.2 1.386 1.000 0 1.000
Nereis.pulsatoria Nototropis.guttatus Oligochaeta
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0 1.000 4.485 0.744 5.602
glms.lvm.zostera.2 3.251 0.941 0.116 1.000 9.386
Paradoneis.harpagonea Parthenina.interstincta
Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 0.616 0 1.000 0.919
glms.lvm.zostera.2 0.179 1.386 1.000 1.662
Parvicardium.exiguum Perinereis.cultrifera
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 0.824 1.000 0.973 1.000
glms.lvm.zostera.2 0.995 0.063 1.000 2.035 0.994
Perioculodes.longimanus Phoronida Phyllodoce.sp.
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 2.111 0.995 1.284 1.000 0
glms.lvm.zostera.2 0 1.000 0.11 1.000 0
Platyhelminthes Platynereis.dumerilii
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.372 1.000 1.257 1.000
glms.lvm.zostera.2 1.000 2.602 0.985 0.374 1.000
Polititapes.aureus Polychaeta.larvae
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.751 1.000 0 1.000
glms.lvm.zostera.2 6.933 0.406 0 1.000
Polydora.ciliata Polygordius.neapolitanus
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 13.8 0.016 0 1.000
glms.lvm.zostera.2 12.416 0.052 0 1.000
Prionospio.cirrifera Protodorvillea.kefersteini
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 28.158 0.002 7.725 0.295
glms.lvm.zostera.2 4.922 0.768 4.879 0.771
Pseudocuma.longicorne Rissoa.membranacea
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.236 1.000 1.543 0.999
glms.lvm.zostera.2 1.386 1.000 1.019 1.000
Rissoa.splendida Salvatoria.clavata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 5.569 0.616 15.955 0.007
glms.lvm.zostera.2 0 1.000 0 1.000
Schistomeringos.rudolphi Sphaerosyllis.hystrix
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 2.389 0.972 0.023 1.000
glms.lvm.zostera.2 8.799 0.222 4.024 0.898
Spio.filicornis Spisula.subtruncata
Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 6.523 0.442 0 1.000
glms.lvm.zostera.2 0.246 1.000 0 1.000
Stenosoma.capito Syllis.gracilis Syllis.hyalina
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 2.178 0.995 4.755 0.709 2.507
glms.lvm.zostera.2 6.661 0.490 17.831 0.004 0
Tellina.tenuis Thracia.phaseolina
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 0.970 1.622 0.999 0.811 1.000
glms.lvm.zostera.2 1.000 1.386 1.000 0 1.000
Tricolia.pullus Tritia.neritea Tritia.reticulata
Dev Pr(>Dev) Dev Pr(>Dev) Dev
glms.lvm.zostera.1
glms.lvm.zostera.3 5.332 0.635 0.632 1.000 0
glms.lvm.zostera.2 1.38 1.000 3.479 0.937 8.301
Turbellaria Upogebia.pusilla
Pr(>Dev) Dev Pr(>Dev) Dev Pr(>Dev)
glms.lvm.zostera.1
glms.lvm.zostera.3 1.000 1.392 1.000 5.769 0.557
glms.lvm.zostera.2 0.264 1.125 1.000 0 1.000
Arguments:
Test statistics calculated assuming uncorrelated response (for faster computation)
P-value calculated using 999 resampling iterations via PIT-trap resampling (to account for correlation in testing.
Well I don’t know how to interpret that.. Have to check the manual first.